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]]
    even = recursive_fft(x[0::2])
    odd = recursive_fft(x[1::2])
    X = [0j] * N
    for k in range(N // 2):
        t = cmath.exp(-2j * math.pi * k / N) * odd[k]
        X[k] = even[k] + t
        X[k + N // 2] = even[k] - t
    return X

def bit_reverse_indices(N):
    bits = N.bit_length() - 1
    order = []
    for i in range(N):
        rev = 0
        v = i
        for _ in range(bits):
            rev = (rev << 1) | (v & 1)
            v >>= 1
        order.append(rev)
    return order

def iterative_fft(x, count_ops=False):
    N = len(x)
    if N & (N - 1) != 0:
        raise ValueError("length must be a power of 2")
    order = bit_reverse_indices(N)
    X = [x[order[i]] for i in range(N)]
    mults = 0
    adds = 0
    size = 2
    while size <= N:
        half = size // 2
        twiddles = [cmath.exp(-2j * math.pi * j / size) for j in range(half)]
        for start in range(0, N, size):
            for j in range(half):
                w = twiddles[j]
                even_part = X[start + j]
                odd_part = X[start + j + half] * w
                mults += 1
                X[start + j] = even_part + odd_part
                X[start + j + half] = even_part - odd_part
                adds += 2
        size *= 2
    if count_ops:
        return X, mults, adds
    return X

print("বিট-রিভার্সাল ক্রম, N=8:", bit_reverse_indices(8))

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_rec = recursive_fft(x)
X_fft_iter = iterative_fft(x)

print()
print(f"{'k':>3} | {'direct DFT':>10} | {'রিকার্সিভ FFT':>13} | {'ইটারেটিভ FFT':>13} | {'|DFT-iter|':>10} | {'|rec-iter|':>10}")
max_diff_dft_iter = 0.0
max_diff_rec_iter = 0.0
for k in range(N):
    d1 = abs(X_dft[k] - X_fft_iter[k])
    d2 = abs(X_fft_rec[k] - X_fft_iter[k])
    max_diff_dft_iter = max(max_diff_dft_iter, d1)
    max_diff_rec_iter = max(max_diff_rec_iter, d2)
    print(f"{k:>3} | {abs(X_dft[k]):>10.6f} | {abs(X_fft_rec[k]):>13.6f} | {abs(X_fft_iter[k]):>13.6f} | {d1:>10.2e} | {d2:>10.2e}")

print()
print(f"সর্বোচ্চ |DFT - ইটারেটিভ FFT|: {max_diff_dft_iter:.3e}")
print(f"সর্বোচ্চ |রিকার্সিভ FFT - ইটারেটিভ FFT|: {max_diff_rec_iter:.3e}")

def counted_dft(x):
    N = len(x)
    mults = 0
    adds = 0
    for k in range(N):
        for n in range(N):
            _ = x[n] * cmath.exp(-2j * math.pi * k * n / N)
            mults += 1
            if n > 0:
                adds += 1
    return mults, adds

print("--- অপারেশন কাউন্ট: direct DFT বনাম FFT ---")
print(f"{'N':>5} | {'DFT মোট':>9} | {'FFT গুণন':>9} | {'FFT যোগ':>8} | {'FFT মোট':>8} | {'স্পিডআপ':>8}")

for N in [8, 16, 32, 64, 128]:
    xN = [math.sin(2 * math.pi * 3 * n / N) for n in range(N)]
    dft_mults, dft_adds = counted_dft(xN)
    dft_total = dft_mults + dft_adds
    _, fft_mults, fft_adds = iterative_fft(xN, count_ops=True)
    fft_total = fft_mults + fft_adds
    speedup = dft_total / fft_total
    print(f"{N:>5} | {dft_total:>9} | {fft_mults:>9} | {fft_adds:>8} | {fft_total:>8} | {speedup:>7.2f}x")