ক্যাপস্টোন — একটি সম্পূর্ণ ডিজিটাল ফিল্টার সিস্টেম এন্ড-টু-এন্ড ডিজাইন ও বিশ্লেষণ

# ========== ধাপ ১: নয়েজ-সহ টেস্ট সিগন্যাল সেটআপ (হাইপোথেটিক্যাল ভাইব্রেশন সেন্সর) ==========
import math, cmath, random

FS = 256.0        # স্যাম্পল রেট, Hz (L56-এর ডিজাইন)
N = 128           # স্যাম্পল সংখ্যা
F_DESIRED = 40.0  # কাঙ্ক্ষিত ভাইব্রেশন টোন, Hz
F_INTERF = 100.0  # ইন্টারফেয়ারেন্স টোন, Hz
A_DESIRED = 1.0
A_INTERF = 0.6
NOISE_SD = 0.3

random.seed(42)
clean = [A_DESIRED * math.sin(2 * math.pi * F_DESIRED * n / FS) for n in range(N)]
interference = [A_INTERF * math.sin(2 * math.pi * F_INTERF * n / FS + 0.7) for n in range(N)]
noise = [random.gauss(0.0, NOISE_SD) for _ in range(N)]
noisy = [clean[n] + interference[n] + noise[n] for n in range(N)]

print(f"fs={FS} Hz, N={N} স্যাম্পল, সময়কাল={N/FS:.3f} সেকেন্ড, বিন-রেজোলিউশন={FS/N:.2f} Hz")
print(f"কাঙ্ক্ষিত: {A_DESIRED}*sin(2*pi*{F_DESIRED}*t), ইন্টারফেয়ারেন্স: {A_INTERF}*sin(2*pi*{F_INTERF}*t+0.7)")
print(f"নয়েজ: gauss(0,{NOISE_SD}), seed=42\n")
for n in range(5):
    print(f" n={n:3d}  clean={clean[n]:8.4f}  interference={interference[n]:8.4f}  noise={noise[n]:8.4f}  noisy={noisy[n]:8.4f}")

# ========== ধাপ ২: হাতে-লেখা DFT দিয়ে স্পেকট্রাল সনাক্তকরণ ==========
def dft(x):
    N = len(x)
    X = []
    for k in range(N):
        s = 0j
        for n in range(N):
            s += x[n] * cmath.exp(-2j * math.pi * k * n / N)
        X.append(s)
    return X

X_noisy = dft(noisy)
mags_noisy = [abs(v) for v in X_noisy]
bin_desired = round(F_DESIRED * N / FS)
bin_interf = round(F_INTERF * N / FS)

print(f"bin_desired (40Hz) = {bin_desired}, bin_interf (100Hz) = {bin_interf}\n")
print(f"{'bin':>4} | {'freq(Hz)':>8} | {'|X|':>10}")
top = sorted(range(1, N // 2), key=lambda k: -mags_noisy[k])[:6]
for k in sorted(top):
    print(f"{k:4d} | {k*FS/N:8.1f} | {mags_noisy[k]:10.4f}")

print(f"\n|X[{bin_desired}]| (40Hz, ফিল্টারের আগে) = {mags_noisy[bin_desired]:.4f}")
print(f"|X[{bin_interf}]| (100Hz, ফিল্টারের আগে) = {mags_noisy[bin_interf]:.4f}")

# ========== ধাপ ৩: FIR ফিল্টার ডিজাইন (L56-এর টুলকিট) + প্রয়োগ ==========
def sinc(x):
    if x == 0.0:
        return 1.0
    return math.sin(math.pi * x) / (math.pi * x)

def design_fir_lowpass(fc, fs, M):
    taps = []
    alpha = (M - 1) / 2
    for n in range(M):
        h_ideal = (2 * fc / fs) * sinc((2 * fc / fs) * (n - alpha))
        w = 0.54 - 0.46 * math.cos(2 * math.pi * n / (M - 1))  # Hamming window
        taps.append(h_ideal * w)
    s = sum(taps)
    return [t / s for t in taps]

FC = 60.0
M_TAPS = 25
h = design_fir_lowpass(FC, FS, M_TAPS)
DELAY = (M_TAPS - 1) // 2

def apply_fir(x, h):
    N = len(x)
    M = len(h)
    y = []
    for n in range(N):
        acc = 0.0
        for k in range(M):
            if n - k >= 0:
                acc += h[k] * x[n - k]
        y.append(acc)
    return y

filtered = apply_fir(noisy, h)

print(f"FIR লো-পাস: fc={FC} Hz, M={M_TAPS} ট্যাপ (Hamming window), group delay={DELAY} স্যাম্পল")
print("প্রথম ৫টি ট্যাপ:", [round(v, 5) for v in h[:5]])
print(f"ট্যাপের সমষ্টি (DC গেইন চেক) = {sum(h):.6f}\n")
for n in range(10, 16):
    print(f" n={n:3d}  noisy={noisy[n]:8.4f}  filtered={filtered[n]:8.4f}")

# ========== ধাপ ৪: আগে/পরে স্পেকট্রাল তুলনা ==========
X_filtered = dft(filtered)
mags_filtered = [abs(v) for v in X_filtered]

print(f"{'freq':>6} | {'|X| আগে':>11} | {'|X| পরে':>10} | {'পরিবর্তন':>8}")
print(f"{F_DESIRED:6.1f} | {mags_noisy[bin_desired]:11.4f} | {mags_filtered[bin_desired]:10.4f} | "
      f"{(mags_filtered[bin_desired]/mags_noisy[bin_desired]-1)*100:7.1f}%")
print(f"{F_INTERF:6.1f} | {mags_noisy[bin_interf]:11.4f} | {mags_filtered[bin_interf]:10.4f} | "
      f"{(mags_filtered[bin_interf]/mags_noisy[bin_interf]-1)*100:7.1f}%")

atten_ratio = mags_noisy[bin_interf] / mags_filtered[bin_interf]
preserve_ratio = mags_filtered[bin_desired] / mags_noisy[bin_desired]
print(f"\nইন্টারফেয়ারেন্স অ্যাটেনুয়েশন ফ্যাক্টর = {atten_ratio:.2f}x")
print(f"কাঙ্ক্ষিত-টোন প্রিজার্ভেশন রেশিও = {preserve_ratio:.4f} ({preserve_ratio*100:.1f}%)")

# ========== ধাপ ৫: SNR, MSE ও correlation মেট্রিক্স (group-delay সংশোধন সহ) ==========
def signal_power(x):
    return sum(v * v for v in x) / len(x)

# SNR (আগে): কাঙ্ক্ষিত সিগন্যাল বনাম (ইন্টারফেয়ারেন্স + নয়েজ)
noise_component_before = [interference[n] + noise[n] for n in range(N)]
P_signal = signal_power(clean)
P_noise_before = signal_power(noise_component_before)
SNR_before_db = 10 * math.log10(P_signal / P_noise_before)

# group-delay সংশোধন করে ফিল্টার করা সিগন্যাল ক্লিন রেফারেন্সের সাথে অ্যালাইন করা
valid_n = range(DELAY, N)
clean_aligned = [clean[n - DELAY] for n in valid_n]
filtered_aligned = [filtered[n] for n in valid_n]
noisy_aligned = [noisy[n] for n in valid_n]

residual_after = [filtered_aligned[i] - clean_aligned[i] for i in range(len(clean_aligned))]
P_signal_aligned = signal_power(clean_aligned)
P_noise_after = signal_power(residual_after)
SNR_after_db = 10 * math.log10(P_signal_aligned / P_noise_after)
SNR_improvement = SNR_after_db - SNR_before_db

mse_before = sum((noisy_aligned[i] - clean_aligned[i]) ** 2 for i in range(len(clean_aligned))) / len(clean_aligned)
mse_after = sum((filtered_aligned[i] - clean_aligned[i]) ** 2 for i in range(len(clean_aligned))) / len(clean_aligned)

def correlation(a, b):
    ma = sum(a) / len(a)
    mb = sum(b) / len(b)
    num = sum((a[i] - ma) * (b[i] - mb) for i in range(len(a)))
    da = math.sqrt(sum((v - ma) ** 2 for v in a))
    db = math.sqrt(sum((v - mb) ** 2 for v in b))
    return num / (da * db)

corr_before = correlation(noisy_aligned, clean_aligned)
corr_after = correlation(filtered_aligned, clean_aligned)

print(f"SNR (আগে)  = {SNR_before_db:.3f} dB")
print(f"SNR (পরে)  = {SNR_after_db:.3f} dB")
print(f"SNR উন্নতি = {SNR_improvement:.3f} dB\n")
print(f"MSE বনাম clean (আগে) = {mse_before:.5f}")
print(f"MSE বনাম clean (পরে) = {mse_after:.5f}")
print(f"MSE হ্রাস ফ্যাক্টর = {mse_before/mse_after:.2f}x\n")
print(f"correlation বনাম clean (আগে) = {corr_before:.4f}")
print(f"correlation বনাম clean (পরে) = {corr_after:.4f}")

# ========== ধাপ ৬: চূড়ান্ত সারাংশ ও সিদ্ধান্ত -- ধাপ ১-৫-এর প্রকৃত ফলাফল থেকে ==========
print("=" * 62)
print("চূড়ান্ত মূল্যায়ন সারাংশ -- ভাইব্রেশন-সেন্সর সিগন্যাল ক্লিনআপ (হাইপোথেটিক্যাল)")
print("=" * 62)
print(f"১. স্পেকট্রাল সনাক্তকরণ: কাঙ্ক্ষিত বিন {bin_desired} (|X|={mags_noisy[bin_desired]:.2f}), "
      f"ইন্টারফেয়ারেন্স বিন {bin_interf} (|X|={mags_noisy[bin_interf]:.2f})")
print(f"২. আগে/পরে স্পেকট্রাল তুলনা: ইন্টারফেয়ারেন্স {atten_ratio:.1f}x কমেছে, "
      f"কাঙ্ক্ষিত টোন {preserve_ratio*100:.1f}% সংরক্ষিত")
print(f"৩. SNR: {SNR_before_db:.2f} dB -> {SNR_after_db:.2f} dB (উন্নতি {SNR_improvement:.2f} dB)")
print(f"৪. MSE বনাম clean: {mse_before:.4f} -> {mse_after:.4f} ({mse_before/mse_after:.1f}x হ্রাস)")
print(f"৫. correlation বনাম clean: {corr_before:.4f} -> {corr_after:.4f}")

if atten_ratio >= 20.0 and preserve_ratio >= 0.80:
    filter_verdict = "ফিল্টার তার ডিজাইন-লক্ষ্য পূরণ করেছে: ইন্টারফেয়ারেন্স দৃঢ়ভাবে দমন করেছে অথচ কাঙ্ক্ষিত টোন সংরক্ষণ করেছে"
else:
    filter_verdict = "ফিল্টার আংশিকভাবে কাজ করেছে, তবে ডিজাইন পুনর্বিবেচনা প্রয়োজন"

if SNR_improvement >= 3.0 and corr_after >= 0.90:
    quality_verdict = "সামগ্রিক সিগন্যাল-মান উল্লেখযোগ্যভাবে উন্নত হয়েছে (SNR ও correlation উভয়েই শক্তিশালী উন্নতি)"
else:
    quality_verdict = "সামগ্রিক সিগন্যাল-মানের উন্নতি সীমিত -- আরও ডিজাইন যাচাই প্রয়োজন"

print(f"\nফিল্টার-কার্যকারিতা সিদ্ধান্ত: {filter_verdict}")
print(f"সামগ্রিক গুণমান সিদ্ধান্ত: {quality_verdict}")
print(f"\nচূড়ান্ত সিদ্ধান্ত: L56-এ ডিজাইন করা পাইপলাইন এই হাইপোথেটিক্যাল সিনারিওতে "
      f"{'কার্যকর প্রমাণিত হয়েছে' if (atten_ratio >= 20.0 and preserve_ratio >= 0.80 and SNR_improvement >= 3.0) else 'আংশিকভাবে কার্যকর হয়েছে'}")