تبدیل فوریه گسسته و نیاز به توابع پنجره
تبدیل فوریه گسسته (Discrete-Time Fourier) یا همان DFT، برای محاسبه طیف فرکانسی یک سیگنال زمان-گسسته مورد استفاده قرار میگیرد. نسخه کارآمد آن به لحاظ محاسباتی، تبدیل فوریه سریع (Fast Fourier Transform ) یا FFT نام دارد که معمولا برای محاسبه DFT استفاده میشود. اما همانطور که خیلی از افراد ممکن است با این وضعیت مواجه شده باشند، استفاده از FFT زمانیکه به تنهایی صورت گیرد معمولا شکل طیف دقیقی را به ما ارائه نمیدهد. علت، پدیدهای است تحت عنوان نشت طیفی (spectral leakage).
نشت طیفی با استفاده از تابع پنجره به همراه تبدیل فوریه گسسته به طور چشمگیری میتواند کاهش یابد. در این مقاله، میخواهیم علت وقوع نشت طیفی را بررسی کرده و سپس توابع پنجرهای را معرفی کنیم که میتواند شکل طیف بدست آمده را بهبود بخشند.
برای یک دنباله زمان-گسسته ، تبدیل فوریه گسسته به شکل زیر تعریف میشود:
که در آن
: طیف فرکانسی گسسته دنباله زمان-گسسته است
: تعداد نمونههای و است
: اندیس زمانی است که معمولا از تا در نظر گرفته میشود.
: اندیس فرکانسی است با همان محدوده تغییرات
از فرمول تعریف DFT مشخص است که این تبدیل به سیگنال با طول محدود اعمال میشود. برای زمان نمونهبرداری ، متغیر زمان گسسته عبارت است از:
برای فرکانس نمونهبرداری ، متغیر فرکانس گسسته عبارت است از:
به طور کلی، مختلط است. برای سیگنال حقیقی، بخش حقیقی نسبت به فرکانس زوج است و بخش موهومی، فرد.
تبدیل DFT میتواند طیف هر سیگنال با طول نمونه یا کمتر را محاسبه کند. از آنجاییکه بسیاری از سیگنالهایی که مورد توجه هستند، نقاط شروع و پایان مشخص و تعریف شدهای ندارند، ما در نهایت باید بخشی از سیگنال با طول را برای محاسبه DFT مورد استفاده قرار دهیم. به عنوان مثال، سیگنال موج سینوسی با فرکانس 8 یا 9 هرتز را در نظر بگیرید، که با فرکانس نمونهبرداری شده است. مقدار را برابر با 64 نمونه فرض میکنیم که از به جای برای نمایش بهتر در مراحل بعدی استفاده کردهایم. کد متلب برای تولید یک موج سینوسی 8 هرتزی و محاسبه DFT آن در پایین آمده است که در آن از تابع fft متلب برای محاسبه DFT طبق معادله مطرح شده در بخشهای قبل استفاده کردهایم. این تابع مقدار طیف گسسته را در نقطه محاسبه میکند.
fs= 128; % Hz sample frequency
Ts= 1/fs; % s sample time
L= 64; % total number of samples
f0= 8; % Hz sine frequency
n= 0:L-1; % time index
x= sin(2*pi*f0*n*Ts)*2/L; % discrete-time sinusoid
X = fft(x,L); % discrete spectrum via DFT
Xmag= abs(X); % magnitude of X
k= 0:L-1; % frequency index
f= k*fs/L; % Hz discrete frequency
به منظور بدست آوردن حداکثر اندازه دامنه طیفی برابر با یک، اندازه دامنه سینوسی با مقدار تغییر مقیاس داده شده است. دقت کنید که برای بردار مختلط ، تابع متلب abs اندازه هر عضو را محاسبه میکند.
موج سینوسی زمان-گسسته در بخش بالایی شکل 1 رسم شده است که در آن واضح است که دقیقا 4 سیکل موج سینوسی در قالب 64 نمونه مشاهده میشود. مقادیر DFT در بخش میانی برای فرکانسهای گسسته از تا رسم شده است. شکل پایینی نیز اندازه را نشان میدهد که بر روی بازه تا رسم شده است. همانطور که رابطه فرکانس گسسته نشان میدهد، فرکانسهای گسسته یا زیربخشهای فرکانسی با فاصله هرتز نسبت به هم کنار یکدیگر قرار گرفتهاند. در شکل 1، جزء فرکانسی تنها در 8 هرتز دیده میشود.
تا کنون به نظر همه چیز خوب است. اما اجازه دهید که با همان کد متلب، طیف فرکانسی موج سینوسی دیگری را محاسبه کنیم:
f0= 9; % Hz
موج سینوسی در بخش بالایی شکل 2 رسم شده است که در اینجا 4.5 سیکل در 64 نمونه دیده میشود. اندازه در بخش پایینی رسم شده است. طیف فرکانسی به نظر ناخوشایند میآید که بر روی زیربخشهای فرکانسی زیادی گسترده شده است: این همان پدیده نشت فرکانسی است که قبلا به آن اشاره شد. اگر به محور فرکانسی دقت کنید، یک مشکل واضح دیده میشود: هیچ زیربخش فرکانسی در 9 هرتز دیده نمیشود. همانطور که خواهیم دید، هیچ راه حل بدون نقص و کاملی نیز برای این مساله وجود ندارد اما میتوان شرایط را بهبود بخشید.
شکل 1- DFT موج سینوسی.
بالا: موج سینوسی زمان-گسسته. وسط: برای تا . پایین: برای تا .
شکل 2- DFT موج سینوسی.
بالا: موج سینوسی زمان-گسسته. وسط: برای تا . پایین: برای تا .
یک پنجره نه چندان مناسب
شکل 3 نمایی از برداشت نمونه از یک موج سینوسی را نشان میدهد. بخش بالایی شکل نمونههای موج سینوسی را نشان میدهد. بخش میانی تابع مستطیلی که 64 نمونه با مقدار واحد دارد را ارائه میکند. بخش پایینی نتیجه ضرب موج سینوسی در تابع پنجره را به صورت نمونه به نمونه و متناظر با هم نشان میدهد.
واضح است که یک پنجره برای ثبت نمونه از موج سینوسی باز میشود و این پنجره تنها یک تابع مستطیلی است. پنجرهگذاری زمانی رخ میدهد که بازه زمانی سیگنال ثبت شده فراتر از رود که در بسیاری از مواقع در سیگنالهای واقعی با آن روبرو میشویم. به طور ریاضی، پنجرهگذاری، یک ضرب عضو به عضو بین موج سینوسی و تابع پنجره است که به نوعی مدولاسیون موج سینوسی با تابع پنجره نیز میتواند تلقی شود.
شکل 3- پنجرهگذاری بوسیله تابع پنجره مستطیلی
بالا: موج سینوسی ثبت شده وسط: تابع پنجره مستطیلی. پایین: موج سینوسی پنجرهگذاری شده
زمانیکه از DFT استفاده میشود، میتوان اثر پنجره مستطیلی را با اضافه کردن صفر به انتهای سیگنال، واضحتر مشاهده کرد (نمونههای با مقدار صفر برای n بیشتر از 63 در شکل 3). این کار باعث میشود که گام فرکانسی گسسته () کوچکتری داشته باشیم که در اینجا تعداد کل نمونهها با احتساب صفرها است.
کد ملتبی که در ادامه آمده است، 64 نمونه موج سینوسی را همچون قبل با فرکانس تولید میکند. سپس صفرها برای رساندن طول سیگنال به به انتهای آن اضافه میشوند. پس از آن DFT این سیگنال جدید محاسبه میشود. (دقت کنید که اگرچه در اینجا به سیگنال صفر اضافه کردهایم اما در صورت استفاده از تابع fft متلب برای سیگنال x با طول کمتر از ، خود تابع به صورت خودکار اضافه کردن صفر را انجام میدهد).
fs= 128; % Hz sample frequency
Ts= 1/fs; % s sample time
L= 64; % number of sinewave samples
f0= 8; % Hz sine frequency
n= 0:L-1; % time index
u= sin(2*pi*f0*n*Ts)*2/L; % discrete-time sinusoid
R= 8; % zero-padding factor
N= R*L; % total samples with zero-padding
x= [u zeros(1,N-L)]; % zero-padded signal
X = fft(x,N); % discrete spectrum via DFT
Xmag= abs(X); % magnitude of X
k= 0:N-1; % frequency index
f= k*fs/N; % Hz discrete frequency
سیگنال با صفرهای اضافه شده در بخش بالایی شکل 4 رسم شده است و اندازه دامنه طیف در بخش پایینی آن، که محور فرکانسی به فرکانسهای 0 تا 32 هرتز () محدود شده است. اجزاء طیفی که به دلیل مدولاسیون موج سینوسی توسط پنجره مستطیلی که در شکل 1 قابل مشاهده نبود، اینجا کاملا مشهود است. (اجزاء DFT با نقطه در شکل 1 توسط دایرههای آبی مشخص شده است). بخش طیفی نزدیک به فرکانس 8 هرتز توسط مدولاسیون گسترده شده است و لوبهای فرعی اطراف آن دیده میشود. همچنین مشاهده میکنیم که در این حالت خاص، طیف سینوسی دقیقا روی یک زیربخش فرکانسی قرار میگیرد و اجزاء طیفی در سایر زیربخشهای فرکانسی DFT با نقطه، صفر هستند.
حال فرض کنید که:
f0= 9; % Hz
سیگنال با صفرهای اضافه شده در بخش بالایی شکل 5 رسم شده است و اندازه دامنه طیف در بخش پایینی آن. حال میبینیم که اجزاء نشت طیفی شکل 2 همگی بر روی قلههای لوبهای فرعی طیف قرار میگیرند. به علاوه، اندازه لوبهای فرعی Xmag خیلی آرام با فرکانس کاهش مییابد. به طور اساسی، پنجره مستطیلی یک درهم ریختگی و بینظمی طیفی ایجاد کرده است.
شکل 4- پنجرهگذاری مستطیلی موج سینوسی با قرار دادن صفر در سیگنال.
بالا: موج سینوسی پنجرهگذاری شده با قرار دادن صفر
پایین: اندازه دامنه طیف موج سینوسی پنجرهگذاری شده
شکل 5- پنجرهگذاری مستطیلی موج سینوسی با قرار دادن صفر در سیگنال.
بالا: موج سینوسی پنجرهگذاری شده با قرار دادن صفر
پایین: اندازه دامنه طیف موج سینوسی پنجرهگذاری شده
طیف پنجره مستطیلی
شکل 4 و 5، طیف فرکانسی موجهای سینوسی پنجرهگذاری شده را نشان داد. حال میخواهیم طیف خود پنجره مستطیلی را بیابیم.
کد متلب زیر یک پنجره با مقدار یک ایجاد کرده و سپس به میزان کافی صفر به انتهای آن اضافه میکند تا طول آن به برسد، به همان شکل که در بخش بالایی شکل 6 نشان داده شده است. اندازه DFT این پنجره مستطیلی صفرگذاری شده تحت عنوان Xmag محاسبه میشود و بعد از آن نیمههای چپ و راست Xmag برای بدست آوردن طیف مرکزی به مرکز صفر با هم جابهجا میشوند. این پاسخ فرکانسی شیفت داده شده با واحد فرکانسی بر حسب در بخش پایینی شکل 6 رسم شده است.
fs= 1; % Hz sample frequency
L= 64; % length of rectangular window
R= 8; % zero-padding factor
N= R*L; % total samples with zero-padding
x= [ones(1,L) zeros(1,N-L)]/L; % rectangular window with zero padding
X = fft(x,N); % discrete spectrum via DFT
Xmag= abs(X); % magnitude of X
k= 0:N-1; % frequency index
f= k*fs/N; % Hz frequency
% spectrum centered at 0 Hz
Xmag_shift = fftshift(Xmag); % swap right half and left half of Xmag
f_shift= f- fs/2; % shift freq range to -fs/2:fs/2
شکل تابع اندازه طیف فرکانسی رسم شده در بخش پایینی شکل 6 به صورت زیر است:
که در آن ، اندیس فرکانسی است، طول پنجره است و طول است که شامل صفرهای اضافه شده نیز میشود. طیف پنجره مستطیلی لوبهای فرعی بزرگی دارد که به آرامی با افزایش فرکانس، کاهش مییابد. همانطور که قبلا دیدیم، این لوبهای فرعی ناخواسته در طیفهای فرکانسی سیگنالهای سینوسی ظاهر میشوند. لوبها اگر در مقیاس لگاریتمی رسم شوند حتی بدتر نیز به نظر میآیند. در نهایت، دقت کنید که پهنای باند لوب اصلی به طور معکوس با L متناسب است: منظور فاصله فرکانسی صفر تا صفر است که برابر است با .
شکل 6- طیف فرکانسی پنجره مستطیلی
بالا: پنجره با صفرهای اضافه شده. پایین: اندازه طیف مرکزی به مرکز 0 هرتز
با مقایسه شکل 4 یا 5 با شکل 6، ممکن است این سوال مطرح شود که چرا با اینکه طیف پنجره متقارن است، طیف فرکانسی موج سینوسی پنجرهگذاری شده چنین نیست. اگر ضرب در حوزه زمان (مدولاسیون) معادل کانولوشن در حوزه فرکانس باشد، آیا نباید کانولوشن پنجره با یک جزء فرکانسی متقارن باشد؟ این موضوع در مورد سیگنالهای زمان-پیوسته صدق میکند ولی برای سیگنال زمان-گسسته قاعده چنین است: ضرب در حوزه زمان معادل کانولوشن دایروی در حوزه فرکانس است. به طور اساسی، اجزاء کانولوشن که در کانولوشن خطی خارج از محدوده 0 تا قرار میگیرند، در نمونه دایروی از انتها بازگشت خورده و به طیف اضافه میشوند که باعث عدم تقارن مشاهده شده در شکلهای 4 و 5 میشود.
طیف بهبود یافته با استفاده از توابع پنجرهای
با توجه به اینکه پنجره مستطیلی چنین آثار مخربی روی نتیجه DFT دارد، این سوال مطرح میشود که: آیا راهی برای بهبود اوضاع وجود دارد؟ راه حل تغییر شکل پنجره یا به عبارتی تابع آن است. ما نیاز به پنجرهای داریم که طیف محدودتری نسبت به پنجره مستطیلی داشته باشد.
بخش بالایی شکل 7، نمونههای یک شکل موج سینوسی را همچون شکل 3 نشان میدهد. بخش میانی، تابع پنجرهای با طول را نشان میدهد که شکل صاف با اندازه دامنه که از تقریبا صفر شروع شده و به همان ختم میشود را دارد. بخش پایینی شکل، نتیجه حاصلضرب موج سینوسی در تابع پنجره را به صورت نمونه به نمونه نشان میدهد. ایده اصلی این بود که پنجره در شکل موج ضرب شود و سپس تبدیل DFT آن محاسبه شود. در ابتدا تاثیر این فرآیند کمی شدید به نظر میرسد زیرا بخش زیادی از نمونههای سیگنال دچار تضعیف میشوند. اما این عمل، در نهایت شدت اثر کمتری نسبت به قطع ناگهانی بخشی از شکل سیگنال توسط پنجره مستطیلی را دارد.
شکل 7- بالا: موج سینوسی ثبت شده. وسط: تابع پنجره. پایین: موج سینوسی پنجرهگذاری شده
تابع پنجره نشان داده شده در شکل 7 که معمولا مورد استفاده قرار میگیرد، هنینگ (Hanning) نام دارد و طبق فرمول زیر تعیین میشود:
به طور مثال برای داریم:
دقت کنید که این تابع توسط تابع متلب hann(L,’periodic’) قابل محاسبه است.
اجازه دهید تا طیف موج سینوسی که توسط پنجره هنینگ پنجرهگذاری شده است را محاسبه کنیم. دوباره، از موج سینوسی با فرکانس 9 هرتز و طول نمونه، با فرکانس نمونهبردای استفاده میکنیم. در کد متلب زیر، یک پنجره هنینگ با طول نمونه ایجاد میشود، سپس موج سینوسی به صورت نمونه به نمونه در تابع پنجره ضرب میشود (این کار با .* در متلب انجام میشود). بعد از آن به موج سینوسی پنجرهگذاری شده ، آنقدر صفر اضافه میکنیم تا بردار با طول 512 نمونه حاصل شود (فاکتور 2 در محاسبه اندازه طیف را تغییر مقیاس میدهد تا مقدار ماکزیمم آن یک شود). موج سینوسی پنجرهگذاری شده که صفر نیز به آن اضافه شده در بخش بالایی شکل 8 نشان داده شده است. این را به خاطر بسپارید که اضافه کردن صفر برای افزایش رزولوشن فرکانسی طیف سیگنال استفاده میشود اما یک پیشنیاز برای محاسبه DFT نیست. در نهایت، تبدیل DFT محاسبه میشود. اندازه طیف حاصل در بخش پایینی شکل 8 آمده است که محور فرکانسی به فرکانسهای 0 تا 32 هرتز () محدود شده است.
fs= 128; % Hz sample frequency
Ts= 1/fs; % s sample time
f0 = 9; % Hz sine frequency
L= 64; % length of window
R= 8; % zero-padding factor
N= R*L; % total samples with zero-padding
n= 0:L-1; % time index prior to zero padding
win= .5*(1-cos(2*pi*n/L)); % hanning window of length L
u= sin(2*pi*f0*n*Ts)*2/L; % discrete-time sine
u_win= win.*u; % windowed sine
x= 2*[u_win zeros(1,N-L)]; % windowed, zero-padded sine
X = fft(x,N); % discrete spectrum via DFT
Xmag= abs(X); % magnitude of X
k= 0:N-1; % frequency index
f= k*fs/N; % Hz frequency
شکل 8- موج سینوسی پنجرهگذاری شده توسط هنینگ با اضافه کردن صفر.
بالا: موج سینوسی پنجرهگذاری شده که به آن صفر اضافه شده است
پایین: اندازه طیف موج سینوسی پنجرهگذاری شده
با مقایسه طیف فرکانسی شکل 8 با نمونه مشابه در شکل 5 ملاحظه میکنیم که لوبهای فرعی خیلی پایینتر قرار میگیرند. هر چند که هنوز طیف نزدیک 9 هرتز به دلیل مدولاسیون توسط تابع پنجره دچار گستردگی شده است (حتی نسبت به پنجره مستطیلی در این حالت عریضتر نیز شده است). متاسفانه این مساله در هر نوع تابع پنجره رخ میدهد.
حال میتوانیم هر دو طیف فرکانسی را در مقیاس دسیبل (dB) با کد زیر رسم کنیم:
XdB= 20*log10(Xmag);
طیف در مقیاس dB در شکل 9 نشان داده شده است. طیف با برچسب “پنجره مستطیلی” همان چیزی است که اگر تنها تبدیل DFT را به موج سینوسی ثبت شده اعمال کنید، دریافت خواهید کرد: این دقیقا طیف موج سینوسی است که توسط تابع پنجره مستطیلی مدوله شده است. طیف با برچسب “پنجره هنینگ” از مدوله کردن موج سینوسی توسط پنجره هنینگ قبل از اعمال تبدیل DFT حاصل میشود. همانطور که قبلا ذکر شد، استفاده از پنجره هنینگ به طور چشمگیری لوبهای فرعی را بهبود بخشید، اما لوب اصلی به طور آزاردهندهای عریض است، با پهنای باند 3dB حدودا 3 هرتز. دقت کنید که پهنای باند لوب اصلی به طور معکوس با L متناسب است. بنابراین، با صرف هزینه به واسطه افزایش میتوانیم این پهنای باند را کاهش دهیم. به طور مثال، شکل 10 طیف موج سینوسی پنجرهگذاری شده برای و را نشان میدهد. حالت با پنجره هنینگ تقریبا شبیه یک طیف سینوسی ایدهآل به نظر میرسد.
یک نکته مهم که تاکنون آموختهایم این است که: فرض کنید که یک سیگنال ناشناخته را ثبت کردهاید، تبدیل DFT نسخه صفر اضافه شده و پنجرهگذاری شده سیگنال توسط تابع هنینگ را نیز محاسبه کردهاید، و طیف فرکانسی در اختیار دارید که مشابه شکل 10 است. ممکن است چنین تصور کنید که لوبهای فرعی دیده شده در شکل جزء شاخصها و ویژگیهای خود سیگنال است. اما حال میدانیم که لوبهای فرعی به دلیل تابع پنجره بوجود آمدهاند و جزء ویژگیهای ذاتی سیگنال نیستند.
طیف شکل 10 از روش اضافه کردن صفر با تبدیل DFT با طول استفاده کرده است. اضافه کردن صفر برای بدست آوردن یک طیف فرکانسی صحیح ضروری نیست: برای مثال، شکل 11 طیف فرکانسی بدون اضافه کردن صفر را نشان میدهد. اینجا، است و گام فرکانسی گسسته است. (بنابراین، انتخاب منجر به قرار گرفتن این فرکانس بین دو زیربخش فرکانسی گسسته میشود). دقت کنید که قلههای طیفهای فرکانسی، ارتفاع کمتری نسبت به قلههای مربوط به طیف حاصل از سیگنال صفر اضافه شده را دارند. این مساله به این دلیل رخ میدهد که در حالت با اضافه کردن صفر روی یک زیربخش فرکانسی مشخص قرار میگیرد ولی برای حالت بدون اضافه کردن صفر چنین اتفاقی نمیافتد. در نهایت، باید به خاطر بسپاریم که دامنههای جانبی پاسخ فرکانسی همان نقاط روی لوبهای فرعی در شکل 10 هستند.
طیف خود توابع پنجره در شکل 12 در مقیاس dB برای حالتهای مستطیلی و هنینگ مقایسه شده است. هر دو تابع دارای پهنای باند متناسب با به طور معکوس هستند. دو ویژگی خیلی مهم تابع پنجره عبارت است از: پهنای باند (پهنای باند نویز) و سطح لوب فرعی. پنجره مستطیلی لوبهای فرعی بلندی دارد که خیلی با سرعت کم با افزایش فرکانس دچار تضعیف میشوند (اولین لوب فرعی در اندازه است). پنجره هنینگ لوبهای فرعی کوتاهتری دارد (اولین لوب فرعی در ) ولی پهنای باند عریضتری دارد. به طور کلی، توابع پنجره با سطح لوب فرعی پایینتر دارای پهنای باند بیشتر هستند. سطح لوب فرعی به طور خاص زمانی مهم است که سعی در نمایش یک سیگنال کوچک در حضور یک سیگنال بزرگ در فرکانس متفاوت داریم.
توابع پنجرهای زیادی با لوبهای فرعی کوتاهتر نسبت به پنجره هنینگ وجود دارند. مثلا پنجره کایزر (Kaiser) این ویژگی را دارد که میتوان سطح لوب فرعی آن را تنظیم کرد.
شکل 9- طیف موجهای سینوسی با استفاده از پنجرههای هنینگ و مستطیلی
شکل 10- طیف موجهای سینوسی با استفاده از پنجرههای هنینگ و مستطیلی
شکل 11- طیف موجهای سینوسی با استفاده از پنجرههای هنینگ و مستطیلی
شکل 12- طیف پنجرههای هنینگ و مستطیلی
منبع: https://www.dsprelated.com/
دیدگاه ها (0)