FIR ফিল্টার — লিনিয়ার ফেজ ও উইন্ডোয়িং দিয়ে ডিজাইন

১ · আদর্শ লো-পাস ফিল্টার ও সিংক ফাংশন

একটি আদর্শ (ideal) লো-পাস ফিল্টারের ফ্রিকোয়েন্সি রেসপন্স একটি "ব্রিক-ওয়াল" আকৃতির — কাটঅফ ফ্রিকোয়েন্সি fc-এর নিচে গেইন ঠিক ১, উপরে ঠিক ০। এই আদর্শ রেসপন্সের ইনভার্স DTFT (M5/L20-এর ধারণা) নিলে যে ইমপালস রেসপন্স পাওয়া যায় তা একটি সুপরিচিত, স্ট্যান্ডার্ড ফলাফল — একটি সিংক (sinc) ফাংশন:

$$h_{ideal}[n] = \frac{\sin(2\pi f_c (n-M_0))}{\pi (n-M_0)}, \qquad h_{ideal}[M_0] = 2f_c$$

এখানে fc normalized কাটঅফ ফ্রিকোয়েন্সি (cycles/sample, অর্থাৎ fc_Hz / fs), আর M0 ফিল্টারের কেন্দ্র (center tap)। n = M0-এ সরাসরি সূত্রটি 0/0 হয়ে যায় বলে সেখানে লিমিট মান 2·fc ব্যবহার করতে হয় — এটি একটি স্ট্যান্ডার্ড, সুপরিচিত গাণিতিক ফলাফল (L'Hôpital-এর নিয়ম দিয়ে প্রমাণযোগ্য)। সমস্যা হলো — এই সিংক ফাংশন উভয় দিকে অসীম পর্যন্ত বিস্তৃত, যা বাস্তবায়ন করা অসম্ভব।

২ · উইন্ডোয়িং পদ্ধতি — একটি বাস্তবায়নযোগ্য ফিল্টার

সমাধান হলো — সিংক ফাংশনকে একটি সসীম দৈর্ঘ্যে কেটে নেওয়া এবং L24-এ শেখা একটি উইন্ডো ফাংশন (এখানে হ্যামিং উইন্ডো) দিয়ে গুণ করা, যা কাটার ফলে হওয়া তীক্ষ্ণ ধার (abrupt edge) মসৃণ করে এবং L24-এর স্পেকট্রাল লিকেজ কমায়:

$$w[n] = 0.54 - 0.46\cos\!\left(\frac{2\pi n}{M-1}\right), \qquad h[n] = h_{ideal}[n]\cdot w[n]$$

চূড়ান্ত ফিল্টার h[n] — n=0 থেকে n=M-1 পর্যন্ত সসীম দৈর্ঘ্যের — এখন একটি বাস্তবায়নযোগ্য FIR ফিল্টার। যেহেতু h_ideal[n] এবং w[n] উভয়ই কেন্দ্র M0=(M-1)/2-এর চারপাশে সিমেট্রিক (h[n] = h[M-1-n]), তাই চূড়ান্ত h[n]-ও সিমেট্রিক — এটাই ঠিক লিনিয়ার ফেজের শর্ত (Type I লিনিয়ার-ফেজ FIR, যখন M বিজোড়)। L32-এ ভিন্ন একটি ডিজাইন পদ্ধতি (ফ্রিকোয়েন্সি-স্যাম্পলিং) দেখা যাবে যেখানে একই সিমেট্রি-শর্ত সরাসরি কাজে লাগে।

৩ · সম্পূর্ণ পাইপলাইন — ডিজাইন করুন, প্রয়োগ করুন, যাচাই করুন

নিচে ধাপে ধাপে সম্পূর্ণ পাইপলাইন চালানো হয়েছে: (ক) fs=200 Hz স্যাম্পল রেটে একটি টেস্ট সিগন্যাল তৈরি — একটি ১০ Hz (নিম্ন) ও একটি ৬০ Hz (উচ্চ) কম্পোনেন্টের যোগফল; (খ) fc=30 Hz কাটঅফ, M=31 ট্যাপের একটি উইন্ডোড-সিংক লো-পাস ফিল্টার ডিজাইন; (গ) কনভোলিউশনের মাধ্যমে ফিল্টার প্রয়োগ (ফিল্টারের গ্রুপ ডিলে — (M-1)/2 = 15 নমুনা — এর জন্য আউটপুট ট্রিম করে ইনপুটের সাথে টাইম-অ্যালাইন করা); (ঘ) একটি ছোট ইনলাইন DFT-ম্যাগনিটিউড ফাংশন দিয়ে ফিল্টারিং-এর আগে ও পরে উভয় ফ্রিকোয়েন্সির ম্যাগনিটিউড পরিমাপ।

import math

fs = 200.0     # স্যাম্পল রেট (Hz)
f_lo = 10.0    # নিম্ন-ফ্রিকোয়েন্সি কম্পোনেন্ট (Hz)
f_hi = 60.0    # উচ্চ-ফ্রিকোয়েন্সি কম্পোনেন্ট (Hz)
N = 60         # টেস্ট সিগন্যালের দৈর্ঘ্য

x = [math.sin(2*math.pi*f_lo*n/fs) + math.sin(2*math.pi*f_hi*n/fs) for n in range(N)]

# ---- (খ) উইন্ডোড-সিংক লো-পাস ফিল্টার ডিজাইন ----
fc = 30.0          # কাটঅফ ফ্রিকোয়েন্সি (Hz)
fc_norm = fc / fs  # normalized কাটঅফ
M = 31             # ট্যাপ সংখ্যা (বিজোড়)
center = (M - 1) // 2

def sinc_lp(d, fc_norm):
    if d == 0:
        return 2 * fc_norm
    return math.sin(2*math.pi*fc_norm*d) / (math.pi*d)

h_ideal = [sinc_lp(n - center, fc_norm) for n in range(M)]

def hamming(M):
    return [0.54 - 0.46*math.cos(2*math.pi*n/(M-1)) for n in range(M)]

w = hamming(M)
h = [h_ideal[n]*w[n] for n in range(M)]

print(f"{'n':>3} | {'h[n]':>10}")
for n in range(M):
    print(f"{n:>3} | {h[n]:>10.6f}")
print(f"\nsum(h) (DC gain) = {sum(h):.6f}")

# ---- (গ) কনভোলিউশনের মাধ্যমে ফিল্টার প্রয়োগ ----
def convolve(a, b):
    ny = len(a) + len(b) - 1
    y = [0.0]*ny
    for n in range(ny):
        s = 0.0
        for k in range(len(b)):
            if 0 <= n-k < len(a):
                s += b[k]*a[n-k]
        y[n] = s
    return y

y_full = convolve(x, h)
y = y_full[center:center+N]   # গ্রুপ ডিলে (center নমুনা) বাদ দিয়ে ইনপুটের সাথে টাইম-অ্যালাইন

# ---- (ঘ) DFT-ম্যাগনিটিউড দিয়ে যাচাই ----
def dft_mag_at(sig, f, fs):
    Ns = len(sig)
    total = sum(sig[n]*complex(math.cos(-2*math.pi*f*n/fs), math.sin(-2*math.pi*f*n/fs)) for n in range(Ns))
    return abs(total) / Ns

mag_in_lo, mag_in_hi = dft_mag_at(x, f_lo, fs), dft_mag_at(x, f_hi, fs)
mag_out_lo, mag_out_hi = dft_mag_at(y, f_lo, fs), dft_mag_at(y, f_hi, fs)

print()
print(f"{'component':>14} | {'|X(f)| আগে':>12} | {'|X(f)| পরে':>12} | {'attenuation':>12}")
for label, mi, mo in [("10 Hz (low)", mag_in_lo, mag_out_lo), ("60 Hz (high)", mag_in_hi, mag_out_hi)]:
    atten_db = 20*math.log10(mo/mi)
    print(f"{label:>14} | {mi:>12.6f} | {mo:>12.6f} | {atten_db:>10.2f} dB")