ডাউনস্যাম্পলিং (ডেসিমেশন)

import math

def dft_magnitude(x, fs):
    """হাতে-লেখা DFT (M5) — ফ্রিকোয়েন্সি বিন k=0..N/2 পর্যন্ত ম্যাগনিচিউড ফেরত দেয়"""
    N = len(x)
    mags = []
    for k in range(N // 2 + 1):
        re = sum(x[n] * math.cos(2 * math.pi * k * n / N) for n in range(N))
        im = -sum(x[n] * math.sin(2 * math.pi * k * n / N) for n in range(N))
        mag = math.sqrt(re**2 + im**2)
        mags.append((k * fs / N, mag))
    return mags

fs = 48.0
N = 48
f_low, f_high = 3.0, 8.0     # দুটি উপাদান: ধীর ৩ Hz, দ্রুত ৮ Hz (অ্যামপ্লিটিউড ০.৮)

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
fs_new = fs / M              # ১২ Hz -- নতুন নাইকুইস্ট সীমা ৬ Hz
x_naive = x[::M]             # নিছক প্রতি চতুর্থ স্যাম্পল রাখা -- কোনো ফিল্টার নেই!

print(f"মূল রেট fs = {fs} Hz, N = {N}  |  ডেসিমেটেড রেট fs' = {fs_new} Hz, N' = {len(x_naive)}")
print(f"নতুন নাইকুইস্ট সীমা: {fs_new/2} Hz  --  ৮ Hz উপাদানটি এই সীমার উপরে!")

# ৮ Hz উপাদানটি নতুন রেটে কোন ফ্রিকোয়েন্সিতে "ফোল্ড" (alias) করবে (L01/L06-এর সূত্র)
k = round(f_high / fs_new)
f_alias = f_high - k * fs_new
print(f"\nপ্রত্যাশিত অ্যালিয়াসড ফ্রিকোয়েন্সি: {f_high} Hz -> {f_alias} Hz (মান হিসেবে {abs(f_alias)} Hz)")

print(f"\n{'freq (Hz)':>10} | {'naive |X[k]|':>13}")
for f, m in dft_magnitude(x_naive, fs_new):
    flag = " -- ভুয়া শক্তি (৮ Hz-এর অ্যালিয়াস)" if abs(f - abs(f_alias)) < 1e-9 and m > 1 else ""
    print(f"{f:>10.2f} | {m:>13.3f}{flag}")

def sinc_lowpass(numtaps, fc):
    """windowed-sinc লো-পাস ফিল্টার (M7/L31-এর পদ্ধতি) -- fc = কাটঅফ, fs-এর ভগ্নাংশ হিসেবে (0..0.5)"""
    M = (numtaps - 1) / 2.0
    h = []
    for n in range(numtaps):
        d = n - M
        if abs(d) < 1e-9:
            hn = 2 * fc                                   # কেন্দ্রে সীমা: sin(x)/x -> 2fc যখন x->0
        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))   # Hamming window
        h.append(hn * w)
    return h

def fir_apply(x, h):
    """ডিসক্রিট কনভোলিউশন (M3/L10) -- ফুল-লেংথ আউটপুট"""
    nx, nh = len(x), len(h)
    y = [0.0] * (nx + nh - 1)
    for n in range(nx):
        for k in range(nh):
            y[n + k] += x[n] * h[k]
    return y

numtaps = 25
fc = (fs_new / 2) / fs        # কাটঅফ = নতুন নাইকুইস্ট (৬ Hz), fs=48-এর ভগ্নাংশ হিসেবে
h = sinc_lowpass(numtaps, fc)

y_full = fir_apply(x, h)
delay = (numtaps - 1) // 2    # লিনিয়ার-ফেজ FIR-এর গ্রুপ ডিলে (M7)
y_aligned = y_full[delay: delay + N]
x_filtered_dec = y_aligned[::M]   # -- এবার ফিল্টার-করা সিগন্যাল থেকে ডেসিমেট

print(f"অ্যান্টি-অ্যালিয়াস ফিল্টার: {numtaps} ট্যাপ, কাটঅফ = {fc*fs} Hz")
print(f"\n{'freq (Hz)':>10} | {'naive |X[k]|':>13} | {'filtered |X[k]|':>16}")
mag_naive = dft_magnitude(x_naive, fs_new)
mag_filt = dft_magnitude(x_filtered_dec, fs_new)
for (f, mn), (_, mf) in zip(mag_naive, mag_filt):
    print(f"{f:>10.2f} | {mn:>13.3f} | {mf:>16.3f}")

reduction_pct = 100 * mag_filt[4][1] / mag_naive[4][1]
preserved_pct = 100 * mag_filt[3][1] / mag_naive[3][1]
print(f"\n৪ Hz ভুয়া শক্তি: naive {mag_naive[4][1]:.3f} -> filtered {mag_filt[4][1]:.3f}  ({reduction_pct:.1f}% অবশিষ্ট)")
print(f"৩ Hz প্রকৃত সিগন্যাল: naive {mag_naive[3][1]:.3f} -> filtered {mag_filt[3][1]:.3f}  ({preserved_pct:.1f}% সংরক্ষিত)")