রেডিক্স-২ FFT অ্যালগরিদম — ডেসিমেশন ইন টাইম

১ · সিগন্যালকে জোড় ও বিজোড় ভাগে ভাঙা

DFT সূত্র থেকে শুরু করা যাক:

$$X[k] = \sum_{n=0}^{N-1} x[n] \, e^{-j2\pi kn/N}$$

এখন যোগফলটিকে জোড় ইনডেক্স (n = 0, 2, 4, ...) আর বিজোড় ইনডেক্স (n = 1, 3, 5, ...) — এই দুই ভাগে ভাগ করা যাক। n = 2m (জোড়) আর n = 2m+1 (বিজোড়) লিখে:

$$X[k] = \underbrace{\sum_{m=0}^{N/2-1} x[2m]\, e^{-j2\pi km/(N/2)}}_{E[k]\ \text{— জোড় স্যাম্পলের } N/2\text{-পয়েন্ট DFT}} \;+\; e^{-j2\pi k/N} \underbrace{\sum_{m=0}^{N/2-1} x[2m+1]\, e^{-j2\pi km/(N/2)}}_{O[k]\ \text{— বিজোড় স্যাম্পলের } N/2\text{-পয়েন্ট DFT}}$$

এখানে E[k] আর O[k] প্রতিটি নিজেই একটি ছোট, স্বাধীন N/2-পয়েন্ট DFT — যথাক্রমে জোড় ও বিজোড় স্যাম্পলের উপর। আর W_N^k = e^{-j2\pi k/N}-কে বলা হয় টুইডল ফ্যাক্টর (twiddle factor)। যেহেতু E[k] ও O[k] পর্যায়ক্রমিক (E[k+N/2] = E[k], একই যুক্তিতে O-এর জন্যও), তাই পুরো N-পয়েন্ট ফলাফল মাত্র দুটো সূত্র থেকেই বের করা যায়:

$$X[k] = E[k] + W_N^k\, O[k], \qquad X[k + N/2] = E[k] - W_N^k\, O[k], \qquad k = 0, \dots, N/2-1$$

এই জোড়া সূত্রকেই বলা হয় বাটারফ্লাই অপারেশন (গ্রাফিক্যালি আঁকলে প্রজাপতির ডানার মতো দেখায় বলে এই নাম)। মূল কথা হলো — N-পয়েন্ট DFT বের করতে দুটো N/2-পয়েন্ট DFT বের করলেই চলে, যদি এই সহজ যোগ-বিয়োগ সূত্র দিয়ে জোড়া লাগানো যায়।

২ · একটি সত্যিকারের রিকার্সিভ ইমপ্লিমেন্টেশন

নিচের কোডে দুটো ফাংশন আছে: direct_dft (M5/L22-এর ডিরেক্ট সূত্র, এই পাঠে যাচাইয়ের জন্য সংক্ষেপে আবার লেখা হলো) আর recursive_fft (উপরের বাটারফ্লাই সূত্র সরাসরি প্রয়োগ করে)। দুটো একই সিগন্যালের উপর চালিয়ে ফলাফল তুলনা করা হচ্ছে:

import cmath, math

def direct_dft(x):
    N = len(x)
    return [sum(x[n] * cmath.exp(-2j * math.pi * k * n / N) for n in range(N)) for k in range(N)]

def recursive_fft(x):
    N = len(x)
    if N == 1:
        return [x[0]]
    if N % 2 != 0:
        raise ValueError("length must be a power of 2")
    even = recursive_fft(x[0::2])
    odd = recursive_fft(x[1::2])
    X = [0j] * N
    for k in range(N // 2):
        twiddle = cmath.exp(-2j * math.pi * k / N) * odd[k]
        X[k] = even[k] + twiddle
        X[k + N // 2] = even[k] - twiddle
    return X

N = 8
fs = 8.0
x = [math.sin(2 * math.pi * 1 * n / fs) + 0.5 * math.cos(2 * math.pi * 2 * n / fs) for n in range(N)]

X_dft = direct_dft(x)
X_fft = recursive_fft(x)

print(f"{'k':>3} | {'direct DFT':>12} | {'recursive FFT':>14} | {'পার্থক্য':>10}")
max_diff = 0.0
for k in range(N):
    d = abs(X_dft[k] - X_fft[k])
    max_diff = max(max_diff, d)
    print(f"{k:>3} | {abs(X_dft[k]):>12.8f} | {abs(X_fft[k]):>14.8f} | {d:>10.2e}")

print()
print(f"সর্বোচ্চ পার্থক্য (৮ বিন মিলিয়ে): {max_diff:.3e}")

N2 = 16
x2 = [math.sin(2 * math.pi * 3 * n / N2) + 0.3 * math.cos(2 * math.pi * 5 * n / N2) for n in range(N2)]
X_dft2 = direct_dft(x2)
X_fft2 = recursive_fft(x2)
max_diff2 = max(abs(a - b) for a, b in zip(X_dft2, X_fft2))
print(f"N=16-এর জন্য সর্বোচ্চ পার্থক্য: {max_diff2:.3e}")

৩ · কেন এই ভাঙন বৈধ — একটি সংক্ষিপ্ত যুক্তি

এই পদ্ধতি কাজ করে কারণ W_N = e^{-j2\pi/N} একটি পর্যায়ক্রমিক (periodic) ও প্রতিসাম্যপূর্ণ সংখ্যা — বিশেষভাবে W_N^{k+N/2} = -W_N^k (কারণ e^{-j\pi} = -1)। এই একটি বীজগাণিতিক সত্যই বাটারফ্লাই সূত্রের + আর - চিহ্নদ্বয়ের পেছনের কারণ — দুটো সূত্রই আসলে একই W_N^k O[k] পদ ব্যবহার করে, শুধু চিহ্ন বদলায়। এই একটি পর্যবেক্ষণই পুরো অ্যালগরিদমের দক্ষতার মূল ভিত্তি: একই গুণফল W_N^k O[k] দুইবার আলাদা করে গণনা করার দরকার নেই।