বায়োমেডিক্যাল সিগন্যাল প্রসেসিং-এ DSP — ECG/EEG পরিচিতি

import math, cmath

fs = 250.0        # স্যাম্পল রেট (Hz)
N = 300           # স্যাম্পল সংখ্যা (সময়কাল ১.২ সেকেন্ড)
dt = 1.0 / fs

def gauss_pulse(t, t0, sigma, amp=1.0):
    return amp * math.exp(-((t - t0) ** 2) / (2 * sigma ** 2))

t0a, t0b = 0.15, 0.75    # দুটি সিন্থেটিক "স্পাইক"-এর কেন্দ্র (সেকেন্ড)
sigma = 0.02              # স্পাইকের প্রস্থ (সেকেন্ড)

# *** সম্পূর্ণ সিন্থেটিক, কাল্পনিক সিগন্যাল -- বাস্তব রোগীর ECG ডেটা নয় ***
clean = [gauss_pulse(n * dt, t0a, sigma) + gauss_pulse(n * dt, t0b, sigma) for n in range(N)]

f_line = 50.0       # পাওয়ার-লাইন ইন্টারফেয়ারেন্স ফ্রিকোয়েন্সি (Hz) -- সিন্থেটিক দূষণ
amp_line = 0.3
noisy = [clean[n] + amp_line * math.sin(2 * math.pi * f_line * 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_clean = dft_magnitude(clean)
mags_noisy = dft_magnitude(noisy)

k50 = round(f_line / (fs / N))
print(f"fs = {fs} Hz, N = {N}, রেজোলিউশন = {fs/N:.4f} Hz/bin, ৫০ Hz-এর বিন k = {k50}")
print()
print(f"{'বিন k':>6} | {'ফ্রিকোয়েন্সি':>12} | {'clean |X[k]|':>13} | {'noisy |X[k]|':>13}")
for k in [2, 4, 6, 8, 10, k50]:
    print(f"{k:>6} | {k*fs/N:>10.2f} Hz | {mags_clean[k]:>13.3f} | {mags_noisy[k]:>13.3f}")

fc = 40.0     # লো-পাস কাট-অফ (Hz) -- ৫০ Hz ইন্টারফেয়ারেন্সের নিচে
M = 20
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]

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(noisy, h)
filtered = filtered_full[M:M + N]     # ফিল্টারের M-স্যাম্পল বিলম্ব বাদ দিয়ে align

mags_filtered = dft_magnitude(filtered)

print(f"ফিল্টার: {taps}-ট্যাপ windowed-sinc lowpass, কাট-অফ = {fc} Hz")
print()
print(f"{'বিন k':>6} | {'ফ্রিকোয়েন্সি':>12} | {'ফিল্টারের আগে':>14} | {'ফিল্টারের পরে':>14}")
for k in [2, 4, 6, 8, 10, k50]:
    print(f"{k:>6} | {k*fs/N:>10.2f} Hz | {mags_noisy[k]:>14.3f} | {mags_filtered[k]:>14.3f}")

reduction = mags_filtered[k50] / mags_noisy[k50]
print()
print(f"৫০ Hz-এ ম্যাগনিটিউড: {mags_noisy[k50]:.3f} -> {mags_filtered[k50]:.3f}  "
      f"(অনুপাত {reduction:.4f}, অর্থাৎ {(1-reduction)*100:.1f}% হ্রাস)")

n_peak = round(t0a * fs)
print(f"প্রথম স্পাইকের কাছাকাছি (n={n_peak}) মান -- clean: {clean[n_peak]:.4f}, "
      f"noisy: {noisy[n_peak]:.4f}, filtered: {filtered[n_peak]:.4f}")