পাওয়ার স্পেকট্রাল ডেনসিটি এস্টিমেশন
এই পাঠে যা শিখবেন
- পাওয়ার স্পেকট্রাল ডেনসিটির সংজ্ঞা ও কেন এটি নয়েজি সিগন্যাল বিশ্লেষণে গুরুত্বপূর্ণ
- পেরিওডোগ্রাম সূত্র ($|X[k]|^2/N$) ও তা M৫-M৬-এর DFT/FFT-এর উপর সরাসরি তৈরি
- Python-এ হাতে-লেখা DFT দিয়ে একটি সত্যিকারের পেরিওডোগ্রাম গণনা ও পিক শনাক্তকরণ
- পেরিওডোগ্রামের সীমাবদ্ধতা (ভ্যারিয়ান্স) এবং Welch's method-এর ধারণাগত সমাধান
১ · পাওয়ার স্পেকট্রাল ডেনসিটি কী
একটি সিগন্যালের পাওয়ার স্পেকট্রাল ডেনসিটি (PSD)Power Spectral Densityএকটি সিগন্যালের গড় পাওয়ার ফ্রিকোয়েন্সি অনুযায়ী কীভাবে বিতরণ হয়ে আছে তা বর্ণনাকারী ফাংশন। বলে দেয় সিগন্যালের মোট পাওয়ার কোন কোন ফ্রিকোয়েন্সিতে "কেন্দ্রীভূত" হয়ে আছে। M২২-এ শেখা DFT সরাসরি একটি নির্দিষ্ট-দৈর্ঘ্যের সিগন্যালের ফ্রিকোয়েন্সি-কম্পোনেন্ট দেয়, কিন্তু PSD-র মূল ব্যবহার হয় যখন সিগন্যালে র্যান্ডম নয়েজ থাকে (M১১-এ বিস্তারিত) — তখন প্রশ্ন হয় "গড়ে কোন ফ্রিকোয়েন্সিতে সত্যিকারের সিগন্যাল লুকিয়ে আছে, নয়েজের ভেতরে?"
সবচেয়ে সাধারণ এস্টিমেটর হলো পেরিওডোগ্রাম:
$$S_{xx}[k] = \frac{|X[k]|^2}{N}$$যেখানে $X[k]$ হলো সিগন্যালের $N$-বিন্দুর DFT (M৫-M৬)। এটি আসলে Parseval-এর উপপাদ্যের একটি সরাসরি প্রয়োগ — সময়-ডোমেইনের এনার্জি ফ্রিকোয়েন্সি-ডোমেইনের এনার্জির সমান, আর $N$ দিয়ে ভাগ করে তা "গড় পাওয়ার প্রতি ফ্রিকোয়েন্সি-বিন" আকারে পাওয়া যায়।
পেরিওডোগ্রাম গণনা করতে নতুন কিছু শিখতে হয় না — শুধু M৫/M৬-এর DFT/FFT নিয়ে তার ম্যাগনিটিউড স্কয়ার করে $N$ দিয়ে ভাগ করলেই হয়।
সত্যিকারের সিগন্যাল ফ্রিকোয়েন্সি বিনগুলোতে পাওয়ার "কেন্দ্রীভূত" থাকে, আর র্যান্ডম নয়েজ প্রায় সমানভাবে সব বিনে ছড়িয়ে থাকে — এই পার্থক্যই পিক শনাক্ত করতে সাহায্য করে।
M২৫-এ শেখা নিয়ম এখানেও প্রযোজ্য — বিন-স্পেসিং হলো $f_s/N$, আর দুটি কাছাকাছি ফ্রিকোয়েন্সি আলাদা করতে প্রকৃত সিগন্যাল-দৈর্ঘ্য (রেকর্ডিং সময়) বাড়াতে হয়, শুধু জিরো-প্যাডিং যথেষ্ট নয়।
২ · পেরিওডোগ্রাম — একটি সত্যিকারের Python ডেমো
নিচের কোডে $N=32$ স্যাম্পলের একটি সিগন্যাল তৈরি করা হয়েছে — $f_s=32$ Hz-এ স্যাম্পল করা একটি $f_0=6$ Hz
সাইন-তরঙ্গ (যা ঠিক bin ৬-এ পড়ে, কারণ বিন-স্পেসিং $f_s/N = 1$ Hz), যার সাথে যোগ করা হয়েছে
random.gauss(0, 0.4) নয়েজ (fixed seed = 7, পুনরুৎপাদনযোগ্যতার জন্য)। M৫-এর সূত্র দিয়ে হাতে-লেখা
DFT ব্যবহার করে পেরিওডোগ্রাম গণনা করা হয়েছে — কোনো NumPy/SciPy নেই।
import math, cmath, random
def dft(x):
N = len(x)
X = []
for k in range(N):
s = 0 + 0j
for n in range(N):
s += x[n] * cmath.exp(-2j * math.pi * k * n / N)
X.append(s)
return X
random.seed(7)
N = 32
fs = 32.0
f0 = 6.0 # প্রকৃত ফ্রিকোয়েন্সি -- ঠিক bin 6-এ পড়ে (fs/N = 1 Hz প্রতি বিন)
noise_std = 0.4
x = [math.sin(2 * math.pi * f0 * n / fs) + random.gauss(0, noise_std) for n in range(N)]
X = dft(x)
psd = [(abs(X[k]) ** 2) / N for k in range(N)] # পেরিওডোগ্রাম: |X[k]|^2 / N
print(f"{'k':>3} | {'freq(Hz)':>9} | {'|X[k]|^2/N':>12}")
for k in range(N // 2 + 1):
freq = k * fs / N
print(f"{k:>3} | {freq:>9.2f} | {psd[k]:>12.4f}")
peak_k = max(range(N // 2 + 1), key=lambda k: psd[k])
noise_bins = [psd[k] for k in range(1, N // 2 + 1) if k != peak_k]
print()
print(f"পিক bin: {peak_k} (ফ্রিকোয়েন্সি {peak_k*fs/N} Hz), পিক মান: {psd[peak_k]:.4f}")
print(f"গড় নয়েজ ফ্লোর (পিক ও DC বাদে): {sum(noise_bins)/len(noise_bins):.4f}")
print(f"সর্বোচ্চ নয়েজ-বিন মান: {max(noise_bins):.4f}")৩ · পেরিওডোগ্রামের সীমাবদ্ধতা — ভ্যারিয়ান্স ও Welch's method
উপরের ডেমোতে নয়েজ-বিনগুলোর মান লক্ষ্য করুন — ০.০০৪৭ থেকে ০.৪৪৪২ পর্যন্ত, অনেক বিস্তৃত। এটাই পেরিওডোগ্রামের পরিচিত দুর্বলতা: এটি PSD-র একটি উচ্চ-ভ্যারিয়ান্স এস্টিমেটর — একই সত্যিকারের PSD থেকে ভিন্ন ভিন্ন নয়েজ-নমুনায় পেরিওডোগ্রাম আউটপুট যথেষ্ট নড়াচড়া করে, এবং $N$ বাড়ালেও এই ভ্যারিয়ান্স কমে না (শুধু ফ্রিকোয়েন্সি-রেজোলিউশন বাড়ে, নির্ভুলতা নয়)।
বাস্তবে এই সমস্যা সমাধানের জন্য Welch's method-এর মতো কৌশল ব্যবহার হয়: সিগন্যালকে ছোট ছোট (সম্ভবত ওভারল্যাপিং) অংশে ভাগ করে প্রতিটি অংশের পেরিওডোগ্রাম আলাদাভাবে গণনা করে, তারপর সবগুলোর গড় নেওয়া হয় — একাধিক স্বতন্ত্র নমুনার গড় নিলে ভ্যারিয়ান্স কমে (ঠিক যেমন পরিসংখ্যানে একাধিক পরিমাপের গড় একক পরিমাপের চেয়ে বেশি নির্ভরযোগ্য)। এই কৌশলটি পরের পাঠের STFT-এর সাথে গঠনগতভাবে সম্পর্কিত — উভয়ই সিগন্যালকে উইন্ডো-করা ছোট অংশে ভাগ করে আলাদাভাবে DFT নেয়, শুধু উদ্দেশ্য ভিন্ন (Welch গড় নেয় একটি একক, আরও নির্ভরযোগ্য স্পেকট্রামের জন্য, STFT প্রতিটি ফ্রেম আলাদা রাখে সময়ের সাথে পরিবর্তন দেখতে)।
পেরিওডোগ্রাম ($|X[k]|^2/N$) হলো PSD এস্টিমেশনের সবচেয়ে সহজ পদ্ধতি, এবং M৫-M৬-এর DFT/FFT-এর উপর সরাসরি তৈরি — উপরের ডেমো দেখাল এটি নয়েজের ভেতরে একটি সত্যিকারের সাইন-তরঙ্গের ফ্রিকোয়েন্সি স্পষ্টভাবে শনাক্ত করতে পারে (পিক নয়েজ ফ্লোরের প্রায় ১৮ গুণ উপরে)। কিন্তু এর উচ্চ ভ্যারিয়ান্সের কারণে বাস্তব অ্যাপ্লিকেশনে Welch's method-এর মতো গড়-করা কৌশল ব্যবহার হয়। এই মডিউলের পরের পাঠে (L৪৫) আমরা DFT-কে সময়ের সাথে সাথে চালিয়ে STFT তৈরি করব — যেখানে সিগন্যাল স্থির নয়, সময়ের সাথে বদলায়।
ভাবনার প্রশ্ন
প্রতিটি প্রশ্ন নিজে কিছুক্ষণ ভাবুন — তারপর "→ উত্তর" চাপুন।
প্র ০১ উপরের ডেমোতে পিক ঠিক bin ৬-এ পড়েছে, কোনো "ফাঁকে" (in-between) নয়। এটি কি কাকতালীয়, নাকি ডিজাইন করা হয়েছিল?
ডিজাইন করা হয়েছিল। বিন-স্পেসিং সবসময় $f_s/N$ — এখানে $32/32 = 1$ Hz প্রতি বিন। যেহেতু সিগন্যালের প্রকৃত ফ্রিকোয়েন্সি $f_0=6$ Hz ঠিক একটি বিনের উপর পড়ে ($6 = 6 \times 1$ Hz), তাই পুরো এনার্জি একটি একক বিনে কেন্দ্রীভূত হয়। যদি $f_0$ দুটি বিনের মাঝামাঝি হতো (যেমন ৬.৫ Hz), তাহলে M২৩/M২৪-এ শেখা স্পেকট্রাল লিকেজের কারণে এনার্জি পার্শ্ববর্তী বিনগুলোতেও ছড়িয়ে পড়ত, আর পিক কম তীক্ষ্ণ দেখাত।
প্র ০২ নয়েজ-বিনগুলোর মান (০.০০৪৭ থেকে ০.৪৪৪২) কেন এত বিস্তৃত, যদিও নয়েজ তাত্ত্বিকভাবে সব ফ্রিকোয়েন্সিতে "সমান শক্তি" রাখার কথা?
random.gauss দিয়ে তৈরি নয়েজের প্রত্যাশিত (expected) PSD সব বিনে সমান, কিন্তু
একটি একক, সসীম-দৈর্ঘ্যের ($N=32$) বাস্তবায়নে এই প্রত্যাশা থেকে যথেষ্ট এলোমেলো বিচ্যুতি হবেই — এটাই ৩ নং
অংশে আলোচিত পেরিওডোগ্রামের উচ্চ ভ্যারিয়ান্স সমস্যা। যদি একই পরীক্ষা বহুবার ভিন্ন
random.seed() দিয়ে চালানো হয় এবং ফলাফলগুলোর গড় নেওয়া হয় (Welch's method-এর মূল ধারণা),
তখন নয়েজ-বিনগুলোর গড় মান একটি সমতল রেখার কাছাকাছি স্থিতিশীল হয়ে আসবে।
প্র ০৩
noise_std-কে ০.৪ থেকে বাড়িয়ে ২.০ করলে পিক শনাক্ত করা কি সহজ হবে না কঠিন?
কঠিন হবে। পিকের মান নির্ভর করে সিগন্যাল-অংশের পাওয়ারের উপর, যা নয়েজ থেকে স্বাধীন, কিন্তু নয়েজ ফ্লোরের
গড় মান নয়েজ-ভ্যারিয়ান্সের সাথে সরাসরি সমানুপাতিক। noise_std বাড়ালে নয়েজ ফ্লোর উঁচু হবে, আর
পিক-থেকে-নয়েজ-ফ্লোর অনুপাত (এটাই মূলত M৪৯-এর SNR ধারণার সাথে সম্পর্কিত) কমে যাবে — চরম ক্ষেত্রে পিকটি
নয়েজের মধ্যে হারিয়ে যেতে পারে।
অনুশীলন
-
চিন্তা করুন: উপরের কোডে
f0 = 6.0-কেf0 = 10.0-এ বদলালে (fs=32, N=32 অপরিবর্তিত), পিক কোন bin-এ পড়বে বলে আপনার ধারণা? বিন-স্পেসিং $f_s/N=1$ Hz ব্যবহার করে হিসেব করুন।বিন-স্পেসিং ১ Hz বলে $f_0=10$ Hz ঠিক bin ১০-এ পড়বে ($10/1 = 10$)। কোড রান করলে bin ১০-এ একটি স্পষ্ট পিক দেখা যাবে, ঠিক যেমন মূল ডেমোতে $f_0=6$ Hz bin ৬-এ পড়েছিল।
-
পরীক্ষা করুন: কোড সেলে
f0-কে ১০.০-এ পরিবর্তন করে Run চেপে আপনার হিসেব যাচাই করুন, এবং লক্ষ্য করুন নতুন পিক মান আগের ৮.১৫৭৬-এর কাছাকাছি কি না।রান করলে bin ১০-এ পিক দেখাবে, আর মান আগের পিকের কাছাকাছি থাকবে (সিগন্যালের অ্যামপ্লিটিউড ও নয়েজ পরিসংখ্যান একই রাখা হয়েছে, শুধু ফ্রিকোয়েন্সি বদলেছে) — তবে
random.gauss-এর কল-সংখ্যা ও ক্রম অপরিবর্তিত থাকায় নয়েজ-নমুনাগুলো একদম অভিন্ন থাকবে, শুধু সিগন্যাল-অংশটি ভিন্ন bin-এ যোগ হবে।
আরও পড়ুন · ABCL TECH-এ আপনার পরবর্তী পদক্ষেপ
- শর্ট-টাইম ফুরিয়ার ট্রান্সফর্ম (STFT) ও স্পেকট্রোগ্রাম পরবর্তী পাঠ যখন সিগন্যালের ফ্রিকোয়েন্সি সময়ের সাথে বদলায়, তখন একক পেরিওডোগ্রাম যথেষ্ট নয় — STFT সিগন্যালকে ছোট ছোট সময়-উইন্ডোতে ভাগ করে প্রতিটির আলাদা স্পেকট্রাম দেখায়।
- নয়েজ ও সিগন্যাল-টু-নয়েজ রেশিও M১১ এই পাঠের পিক-বনাম-নয়েজ-ফ্লোর তুলনাটি M১১-এ SNR-এর একটি সুনির্দিষ্ট, গাণিতিক সংজ্ঞায় পরিণত হবে।
- Numerical Methods কোর্স সহোদর কোর্স সাধারণ সংখ্যাগত অ্যালগরিদম ও ত্রুটি বিশ্লেষণ — এই কোর্স সেই একই হাতে-লেখা দর্শনের উপর সিগন্যাল-প্রসেসিং-নির্দিষ্ট দিকটি যোগ করে।