অডিও প্রসেসিং-এ DSP

import math, cmath

fs = 400.0   # স্যাম্পল রেট (Hz) -- শুধু ডেমো দ্রুত রাখতে ছোট রাখা হয়েছে
N = 64       # স্যাম্পল সংখ্যা
f0, f1, f2 = 50.0, 100.0, 150.0   # ফান্ডামেন্টাল + ২য় + ৩য় হারমোনিক

# *** সম্পূর্ণ সিন্থেটিক, শিক্ষামূলক টোন -- কোনো রেকর্ড করা অডিও নয় ***
tone = [
    1.00 * math.sin(2 * math.pi * f0 * n / fs)
    + 0.50 * math.sin(2 * math.pi * f1 * n / fs)
    + 0.25 * math.sin(2 * math.pi * f2 * n / fs)
    for n in range(N)
]

def dft_magnitude(x):
    n_samples = len(x)
    mags = []
    for k in range(n_samples):
        s = 0 + 0j
        for n in range(n_samples):
            s += x[n] * cmath.exp(-2j * math.pi * k * n / n_samples)
        mags.append(abs(s))
    return mags

mags = dft_magnitude(tone)

print("প্রথম ৮টি স্যাম্পল:", [round(v, 4) for v in tone[:8]])
print()
print(f"fs = {fs} Hz, N = {N} স্যাম্পল, ফ্রিকোয়েন্সি রেজোলিউশন = {fs/N} Hz/bin")
print()
print(f"{'k':>3} | {'ফ্রিকোয়েন্সি (Hz)':>18} | {'|X[k]|':>10}")
for k in range(N // 2 + 1):
    if mags[k] > 1.0:
        print(f"{k:>3} | {k*fs/N:>18.2f} | {mags[k]:>10.3f}")

# windowed-sinc লো-পাস ফিল্টার -- কাট-অফ ১২০ Hz (M7-এর পদ্ধতি)
fc = 120.0
M = 16                 # ফিল্টার অর্ধ-দৈর্ঘ্য -> মোট ট্যাপ = 2M+1
taps = 2 * M + 1
fc_norm = fc / fs

h = []
for i in range(-M, M + 1):
    if i == 0:
        hv = 2 * fc_norm
    else:
        hv = math.sin(2 * math.pi * fc_norm * i) / (math.pi * i)
    w = 0.54 - 0.46 * math.cos(2 * math.pi * (i + M) / (taps - 1))   # Hamming window
    h.append(hv * w)

dc_gain = sum(h)
h = [c / dc_gain for c in h]     # DC গেইন ১-এ নরমালাইজ

def convolve(x, h):
    y = [0.0] * (len(x) + len(h) - 1)
    for n in range(len(x)):
        for k in range(len(h)):
            y[n + k] += x[n] * h[k]
    return y

filtered_full = convolve(tone, h)
filtered = filtered_full[M:M + N]     # ফিল্টারের M-স্যাম্পল বিলম্ব বাদ দিয়ে align

mags_after = dft_magnitude(filtered)

print(f"ফিল্টার: {taps}-ট্যাপ windowed-sinc lowpass, কাট-অফ = {fc} Hz")
print()
print(f"{'উপাদান':>13} | {'ফ্রিকোয়েন্সি':>13} | {'আগে |X[k]|':>12} | {'পরে |X[k]|':>12} | {'অনুপাত':>8}")
labels = ["ফান্ডামেন্টাল", "২য় হারমোনিক", "৩য় হারমোনিক"]
for label, k in zip(labels, [8, 16, 24]):
    ratio = mags_after[k] / mags[k]
    print(f"{label:>13} | {k*fs/N:>10.1f} Hz | {mags[k]:>12.3f} | {mags_after[k]:>12.3f} | {ratio:>8.4f}")