পলিফেজ ফিল্টার ও এফিসিয়েন্ট মাল্টিরেট স্ট্রাকচার

import math

def sinc_lowpass(numtaps, fc):
    M = (numtaps - 1) / 2.0
    h = []
    for n in range(numtaps):
        d = n - M
        if abs(d) < 1e-9:
            hn = 2 * fc
        else:
            hn = math.sin(2 * math.pi * fc * d) / (math.pi * d)
        w = 0.54 - 0.46 * math.cos(2 * math.pi * n / (numtaps - 1))
        h.append(hn * w)
    return h

def convolve_counted(x, h):
    """ডিসক্রিট কনভোলিউশন (M3/L10) -- সাথে প্রকৃত গুণন-অপারেশন গুনে রাখা"""
    nx, nh = len(x), len(h)
    y = [0.0] * (nx + nh - 1)
    mults = 0
    for n in range(nx):
        for k in range(nh):
            y[n + k] += x[n] * h[k]
            mults += 1
    return y, mults

def x_get(x, j):
    """রেঞ্জের বাইরের ইনডেক্সকে শূন্য ধরা (zero-padding), L40-এর মতোই"""
    return x[j] if 0 <= j < len(x) else 0.0

# -- L40-এর হুবহু একই টেস্ট সিগন্যাল ও ফিল্টার --
fs, N = 48.0, 48
f_low, f_high = 3.0, 8.0
x = [math.sin(2*math.pi*f_low*n/fs) + 0.8*math.sin(2*math.pi*f_high*n/fs) for n in range(N)]

M = 4
numtaps = 25
fc = (fs/M/2) / fs
h = sinc_lowpass(numtaps, fc)
delay = (numtaps - 1) // 2

# ---- পদ্ধতি ১: সরাসরি (L40-এর মতো) -- পুরো কনভোলিউশন, তারপর ডেসিমেট ----
y_full, mults_direct = convolve_counted(x, h)
y_aligned = y_full[delay: delay + N]
y_direct_dec = y_aligned[::M]
n_out = len(y_direct_dec)

# ---- পদ্ধতি ২: পলিফেজ -- M-টি ফেজ-সাব-ফিল্টার, প্রতিটি নিম্ন রেটে ----
h_phases = [h[p::M] for p in range(M)]
for p in range(M):
    print(f"h_phase[{p}] দৈর্ঘ্য: {len(h_phases[p])} ট্যাপ")

y_poly = [0.0] * n_out
mults_poly = 0
for p in range(M):
    h_p = h_phases[p]
    c = delay - p
    j_min = -(len(h_p) - 1)                                  # সাব-ফিল্টারের জন্য প্রয়োজনীয় "ইতিহাস"
    e_p = [x_get(x, M*j + c) for j in range(j_min, n_out)]    # ফেজ p-এর নিম্ন-রেট ইনপুট সেগমেন্ট
    branch_out, branch_mults = convolve_counted(e_p, h_p)
    mults_poly += branch_mults
    offset = len(h_p) - 1
    for m in range(n_out):
        y_poly[m] += branch_out[m + offset]

print(f"\n{'m':>3} | {'সরাসরি (L40)':>13} | {'পলিফেজ':>13} | {'পার্থক্য':>10}")
max_diff = 0.0
for m in range(n_out):
    diff = abs(y_direct_dec[m] - y_poly[m])
    max_diff = max(max_diff, diff)
    print(f"{m:>3} | {y_direct_dec[m]:>13.6f} | {y_poly[m]:>13.6f} | {diff:>10.2e}")

print(f"\nসর্বোচ্চ পার্থক্য: {max_diff:.2e}  (ফ্লোটিং-পয়েন্ট নির্ভুলতার সীমা -- ব্যবহারিকভাবে অভিন্ন)")
print(f"\nগুণন-অপারেশন -- সরাসরি: {mults_direct},  পলিফেজ: {mults_poly},  অনুপাত: {mults_direct/mults_poly:.2f}x")

# একই ফিল্টার-আকৃতি (h_phases-এর সাব-ফিল্টার দৈর্ঘ্য), বিভিন্ন সিগন্যাল দৈর্ঘ্য N
h_lens = [len(hp) for hp in h_phases]   # [7, 6, 6, 6], উপরের সেল থেকে

print(f"{'N':>6} | {'n_out':>6} | {'mults_direct':>13} | {'mults_poly':>11} | {'অনুপাত':>8}")
for N_test in [48, 96, 480, 4800]:
    n_out_test = N_test // M
    mults_direct_test = N_test * numtaps
    mults_poly_test = sum((n_out_test + hl - 1) * hl for hl in h_lens)
    ratio = mults_direct_test / mults_poly_test
    print(f"{N_test:>6} | {n_out_test:>6} | {mults_direct_test:>13} | {mults_poly_test:>11} | {ratio:>7.2f}x")

print(f"\nতত্ত্বীয় সর্বোচ্চ সাশ্রয় (M): {M}x -- N বাড়ার সাথে সাথে অনুপাত এই সীমার কাছাকাছি পৌঁছায়")