تفسیر تبدیل فوریه سریع (FFT)، تبدیل فوریه گسسته (DFT) مختلط، بازه‌های فرکانسی و fftshift

اغلب با مساله تبدیل یک سیگنال حوزه زمان به حوزه فرکانس و بالعکس مواجه میشویم. تبدیل فوریه یک ابزار عالی برای دستیابی به این تبدیل است و به طور گسترده در بسیاری از کاربردها استفاده میشود. در پردازش سیگنال، یک سیگنال حوزه زمان میتواند پیوسته یا گسسته باشد و همچنین میتواند متناوب یا نامتناوب نیز باشد. این دسته‌بندی باعث ایجاد ۴ نوع تبدیل فوریه میشود.

جدول ۱- انواع چهارگانه تبدیل فوریه

از جدول ۱ میتوان فهمید که زمانیکه سیگنال در یک حوزه گسسته است، در حوزه دیگر متناوب خواهد بود. به طور مشابه، اگر سیگنال در یک حوزه پیوسته باشد، در حوزه دیگر نامتناوب خواهد بود. برای سادگی بیشتر اجازه دهید که به مساله فرمولهای خاص هر کدام از تبدیلهای بالا نپردازیم. بحث خود را به DFT محدود میکنیم که به عنوان بخشی از بسته‌های نرم‌افزاری همچون متلب (MATLAB)، پایتون (python – Scipy) و غیره در دسترس است. البته میتوانیم تبدیلهای دیگر را توسط DFT تقریب بزنیم.

نسخه‌های حقیقی و مختلط تبدیلها

برای هر کدام از تبدیلهای لیست شده در بالا، یک نسخه حقیقی و یک نسخه مختلط وجود دارد. نسخه حقیقی تبدیل، اعداد حقیقی را به عنوان ورودی میگیرد و دو مجموعه از نقاط حوزه فرکانس تحویل میدهد که یک مجموعه، ضرایب روی پایه تابع کسینوسی و مجموعه دیگر، ضرایب روی تابع پایه سینوسی را نمایش میدهد. نسخه مختلط تبدیلها، فرکانسهای مثبت و منفی را در یک آرایه نشان میدهند. نسخه‌های مختلط، انعطاف‌پذیر هستند و میتوانند هر دو نوع سیگنال با مقدار حقیقی و مختلط را پردازش کنند. شکل زیر، اختلاف بین DFT حقیقی و مختلط را به نمایش میگذارد.

شکل ۱- DFT حقیقی و مختلط

DFT‌ حقیقی

DFT‌ به صورت N-نقطه‌ای حقیقی را در نظر بگیرید. N‌ نمونه از شکل موج x[n] در حوزه زمان با مقدار حقیقی را میگیرد و دو آرایه به طول N/2+1‌ برای هر کدام که روی توابع کسینوسی و سینوسی نگاشت شده‌اند به ترتیب تحویل میدهد:

$\begin{aligned} X_{re}[k] = \frac{2}{N} \sum_{n=0}^{N-1} x[n]. cos \left(\frac{2 \pi k n}{N} \right) \\ \\ X_{re}[k] =- \frac{2}{N} \sum_{n=0}^{N-1} x[n]. sin \left(\frac{2 \pi k n}{N} \right) \end{aligend}$

در اینجا، اندیس حوزه زمان، n، در محدوده $0 \to N$ تغییر میکند و اندیس حوزه فرکانس، k، در محدوده $0 \to N/2$ .

سیگنال حوزه زمان با مقدار حقیقی میتنواند از زوج DFT حقیقی به شکل زیر،‌ سنتز شود:

$</span> x[n] = \sum_{k=0}^{N/2} X_{re}[K] cos \left(\frac{2 \pi k n}{N} \right) - X_{im}[K].sin \left(\frac{2 \pi k n}{N} \right)$

اخطار: زمانیکه از معادله سنتز استفاده میشود، مقادیر $X_{re}[0]$ و $X_{re}[N/2]$ باید بر دو تقسیم شوند. این مساله به دلیل آن است که ما تحلیل را به مقادیر حقیقی محدود میکنیم. این نوع مسائل با استفاده از نسخه مختلط DFT قابل اجتناب است.

DFT‌ مختلط

DFT با N-نقطه مختلط را در نظر بگیرید. N نمونه از مقادیر مختلط شکل موج حوزه زمان x[n] را میگیرد و یک آرایه با طول N‌ به شکل زیر تولید میکند:

$X[k] = \sum_{n=0}^{N-1} x[n]e^{-j2 \pi kn/N}$

مقادیر آرایه به شکل زیر تفسیر میشوند:

X[0] جزء فرکانس DC را نشان میدهد
N/2 المان بعدی، اجزاء فرکانس مثبت هستند که X[N/2] فرکانس نایکوئیست است (که برابر با نصف فرکانس نمونه‌برداری است).
N/2-1‌ المان بعدی، اجزاء فرکانس منفی هستند (دقت کنید که اجزاء فرکانس منفی، در فازورهایی هستند که در جهت معکوس میچرخند و به طور دلخواه بر اساس نوع کاربرد قابل حذف هستند).

معادله سنتز متناظر (x[n] را از نمونه‌های حوزه فرکانس X[k] بازیابی میکند) عبارت است از:

$x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k]e^{j2 \pi kn/N}$

از این معادلات مشخص است که DFT حقیقی از طریق نگاشت سیگنال روی توابع سینوسی و کسینوسی محاسبه میشوند. در حالیکه، DFT مختلط، سیگنال ورودی را روی توابع پایه نمایی نگاشت میدهد (فرمول اویلر این مفاهیم را به هم متصل میکند).

زمانیکه سیگنال ورودی در حوزه زمان دارای مقادیر حقیقی است، DFT مختلط در حین محاسبه بخش موهومی را با مقادیر صفر، پر میکند (این موضوع انعطاف‌پذیری تبدیل مختلط را نشان میدهد و از مشکلی که در بالا برای DFT‌ حقیقی اشاره شد، جلوگیری میکند). شکل زیر نشان میدهد که چگونه نتایج خام FFT در متلب که برای محاسبه DFT‌ به کار میرود، تفسیر میشوند. جزییات بیشتر با یک مثال در ادامه مورد بحث قرار خواهد گرفت.

شکل ۲- تفسیر فرکانسها در خروجی DFT مختلط

تبدیل فوریه سریع (FFT)

تابع fft در متلب، الگوریتمی است که در سال 1965 توسط J.W.Cooley و J.W.Tuckey برای محاسبه تبدیل DFT با راندمان بهتر منتشر شد. این روش از ساختار ویژه DFT زمانیکه طول یسگنال توانی از عدد ۲ است، بهره میبرد و در این شرایط پیچیدگی محاسباتی به طور قابل توجهی کاهش می‌یابد. معمولا طول FFT به صورت توانی از ۲ در نظر گرفته میشود که به آن FFT‌ در مبنای ۲ یا 2FFT میگویند که از فاکتورهای پیچیدگی (twiddle factor) استفاده میکند. طول FFT میتواند عدد فرد نیز باشد همانطور که در روش FFT‌ با فاکتورهای اول پیاده‌سازی میشود که در آن طول FFT به دو عدد اول فاکتور میشود.

FFT به طور وسیع در بسته‌های نرم‌افزاری همچون MATLAB, Scipy‌ و غیره در دسترس است. در این نرم‌افزارها، FFT نسخه مختلط تبدیل DFT‌ را پیاده‌سازی میکند. پیاده‌سازی FFT‌ در متلب، DFT‌ مختلطی را محاسبه میکند که خیلی شبیه معادلات بالا است به غیر از فاکتور مقیاس‌دهی. برای مقایسه، پیاده‌سازی FFT‌ در متلب در فرمولهای زیر آمده است:

$\begin{aligned} X[k] = \sum_{n=0}^{N-1} x[n] e^{-j 2 \pi kn/N} \\ X[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j 2 \pi kn/N} \end{aligend}$

دستورهای متلب برای پیاده‌سازی معادلات بالا به ترتیب fft و ifft‌ هستند. نحوه استفاده از این توابع در زیر آمده است:

X = fff(x,N) % compute X[k]

x = ifft(X,N) % compute x[n]

تفسیر نتایج FFT

بیایید فرض کنیم سیگنال x[n] در حوزه زمان، یک سیگنال کسینوسی با فرکانس $f_c=10 Hz$ است که با نرخ فرکانسی $f_s=32*f_c$ برای استفاده در رایانه نمونه‌برداری شده است.

fc=10;%frequency of the carrier

fs=32*fc;%sampling frequency with oversampling factor=32

t=0:1/fs:2-1/fs;%2 seconds duration

x=cos(2*pi*fc*t);%time domain signal (real number)

subplot(3,1,1); plot(t,x);hold on; %plot the signal

title(‘x[n]=cos(2 \pi 10 t)’); xlabel(‘n’); ylabel(‘x[n]’);

شکل ۳- یک بازه ۲ ثانیه‌ای از موج کسینوسی با فرکانس ۱۰ هرتز

حال تبدیل FFT با N=256‌ نقطه را در نظر بگیرید که طول آن در واقع توان ۸ام عدد ۲ است.

N=256; %FFT size

X = fft(x,N);%N-point complex DFT

%output contains DC at index 1, Nyquist frequency at N/2+1 th index

%positive frequencies from index 2 to N/2

%negative frequencies from index N/2+1 to N

نکته:‌طول FFT باید به اندازه‌ای که تمام طول سیگنال ورودی را پوشش دهد، باشد. اگر N کمتر از طول سیگنال ورودی باشد، سیگنال ورودی باید برای محاسبه در FFT‌ کوتاه یا اصلاحاً خلاصه‌سازی (truncate) شود. برای مثال ما، مدت زمان موج کسینوسی مدت زمان ۲ ثانیه است و شامل ۶۴۹ نقطه است (یک موج با فرکانس ۱۰ هرتز که با فاکتور فرانمونه‌برداری ۳۲، نمونه‌برداری شده دارای $2 \times 32 \times 10 = 640$ نمونه در بازه زمانی ۲ ثانیه است. از آنجا که سیگنال ورودی ما متناوب است، میتوانیم با خیال راحت از N=256 برای طول FFT‌ استفاده کنیم، البته در هر صورت در حین محاسبه FFT‌، طول سیگنال گسترش خواهد یافت.

با توجه به اینکه در متلب،‌ اندیس بردار از شماره ۱ شروع میشود، جزء DC تجزیه FFT در اندیس ۱ ظاهر میشود:

>X(1)

1.1762e-14 (approximately equal to zero)

این موضوع خیلی ساده است. دقت کنید که اندیس مقادیر خام FFT، اعداد صحیح در بازه $1 \to N$ هستند. نیاز به انجام یک پردازش برای تبدیل این اعداد صحیح به فرکانسها داریم. اینجا است که فرکانس نمونه‌برداری وارد عمل میشود.

هر نقطه/بازه در آرایه خروجی FFT، توسط رزولوشن فرکانسی $\delta f$ از هم جدا میشوند که به شکل زیر محاسبه میشود:

$\delta f = \frac{f_s}{N}$

که $f_s$ فرکانس نمونه‌برداری است و $N$ اندازه FFT مورد نظر است. بنابراین، برای مثال ما، هر نقطه از آرایه توسط فاصله رزولوشن فرکانسی زیر از هم جدا شده است:

$\delta f = \frac{f_s}{N} = \frac{32*f_c}{256} = \frac{320}{256} = 1.25 Hz$

حال، سیگنال کسینوسی ۱۰ هرتزی، یک جهش فرکانسی در نمونه ۸ام ( $10/1.25=8$ ) خواهد داشت که در اندیس شماره ۹ در متلب قرار دارد.

> abs(X(8:10)) %display samples 7 to 9

ans = 0.0000 128.0000 0.0000

بنابراین، از رزولوشن فرکانسی میتوان تمام محور فرکانسی را به شکل زیر محاسبه کرد:

%calculate frequency bins with FFT

df=fs/N %frequency resolution

sampleIndex = 0:N-1; %raw index for FFT plot

f=sampleIndex*df; %x-axis index converted to frequencies

حال، میتوانیم اندازه مقادیر FFT را بر حسب فرکانس به شکل زیر رسم کنیم:

subplot(3,1,2); stem(sampleIndex,abs(X)); %sample values on x-axis

title(‘X[k]’); xlabel(‘k’); ylabel(‘|X(k)|’);

subplot(3,1,3); plot(f,abs(X)); %x-axis represent frequencies

title(‘X[k]’); xlabel(‘frequencies (f)’); ylabel(‘|X(k)|’);

شکل زیر، محور فرکانس و اندیس نمونه‌ها را همانطور که در خروجی FFT داریم،‌ نمایش میدهد:

شکل ۴- پاسخ اندازه خروجی FFT که بر حسب اندیس آنها رسم شده است (بالا) و فرکانسهای محاسبه شده متناظر (در پایین).

بعد از اینکه محور فرکانس با توجه به فرکانس نمونه‌برداری به شکل صحیح تبدیل شد، متوجه میشویم که سیگنال کسینوسی یک جهش در فرکانس ۱۰ هرتز دارد. علاوه بر این، یک جهش نیز در نمونه ۲۵۶-۸=۲۴۸ دارد که متعلق به بخش فرکانس منفی است. از آنجا که ما به ماهیت سیگنال آگاه هستیم، به طور دلخواه میتوان فرکانس منفی را در نظر نگرفت. نمونه در فرکانس نایکوئیست ( $f_s/2$ )، مرز بین فرکانسهای منفی و مثبت را مشخص میکند.

> nyquistIndex=N/2+1;

> X(nyquistIndex-2:nyquistIndex+2).’

ans = 1.0e-13 *

-0.2428 + 0.0404i

-0.1897 + 0.0999i

-0.3784ب

-0.1897 – 0.0999i

-0.2428 – 0.0404i

دقت کنید که اعداد مختلط حول اندیس نایکوئیست به صورت مزدوج هستند و به نوعی فرکانسهای به ترتیب مثبت و منفی را نمایش میدهند.

FFTshift

از شکل بالا مشخص است که محور فرکانس از مقدار DC‌ شروع شده و به دنبال آن فرکانسهای مثبت و سپس فرکانسهای منفی می‌آیند. برای ترتیب مناسب بر روی محور x، میتوان از تابع fftshift‌ متلب استفاده کرد که ترتیب فرکانسها را به شکل صحیح درمی‌آورد:‌ فرکانسهای منفی -> فرکانس DC‌ -> فرکانسهای مثبت. تابع fftshift زمانیکه N‌ عدد فرد است باید با دقت لازم استفاده شود.

برای N زوج، ترتیب ارائه شده توسط fftshift به شکل زیر است (دقت کنید که تمامی اندیسها در زیر متناظر با اندیس متلب است):

X[1] جزء فرکانس DC را نشان میدهد
X[2] تا X[N/2] اجزاء فرکانس مثبت هستند.
X[N/2+1] فرکانس نایکوئیست ( $F_s/2$ ) است که بین فرکانسهای مثبت و منفی، مشترک است. این فرکانس را بخشی از فرکانسهای منفی در نظر میگیریم تا با تابع fftshift‌ تناظر درستی داشته باشیم.
X[N/2+1] تا X[N] به عنوان اجزاء‌ فرکانسی منفی در نظر گرفته میشوند.

تابع fftshift، جزء DC را به مرکز طیف انتقال یا شیفت میدهد. این مساله که فرکانس نایکوئیست در اندیس شماره N/2+1 در متلب، بین هر دو بخش فرکانس مثبت و منفی مشترک است، مهم است. دستور fftshift فرکانس نایکوئیست را در بخش فرکانسهای منفی قرار میدهد. شکل زیر این مساله را نشان میدهد:

شکل ۵- نقش fftshift در مرتب کردن فرکانسها

بنابراین، زمانیکه N‌ زوج است، محور فرکانسی مرتب شده به شکل زیر است:

$f = \delta f \left[-\frac{N}{2}:1:\frac{N}{2}-1 \right] = \frac{f_s}{N} \left[-\frac{N}{2}:1:\frac{N}{2}-1 \right]$

زمانیکه N‌ فرد است، محور فرکانسی مرتب شده باید به شکل زیر تنظیم شود:

$f = \delta f \left[-\frac{N-1}{2}:1:\frac{N+1}{2}-1 \right] = \frac{f_s}{N} \left[-\frac{N-1}{2}:1:\frac{N+1}{2}-1 \right]$

کد کوتاه زیر، روش fftshift‌ را از هر دو روش دستی و دستور متلب محاسبه میکند. نتایج از طریق روی هم گذاشتن آنها در یک نمودار، نشان داده شده است. شکل نشان میدهد که هر دو روش در توافق با هم هستند و نتایج یکسانی دارند.

%two-sided FFT with negative frequencies on left and positive frequencies on right

%following works only if N is even, for odd N see equation above

X1 = [(X(N/2+1:N)) X(1:N/2)]; %order frequencies without using fftShift

X2 = fftshift(X);%order frequencies by using fftshift

df=fs/N; %frequency resolution

sampleIndex = -N/2:N/2-1; %raw index for FFT plot

f=sampleIndex*df; %x-axis index converted to frequencies

%plot ordered spectrum using the two methods

figure;subplot(2,1,1);stem(sampleIndex,abs(X1));hold on; stem(sampleIndex,abs(X2),’r’)

%sample index on x-axis

title(‘Frequency Domain’); xlabel(‘k’); ylabel(‘|X(k)|’);%x-axis represent sample index

subplot(2,1,2);stem(f,abs(X1)); stem(f,abs(X2),’r’) %x-axis represent frequencies

xlabel(‘frequencies (f)’); ylabel(‘|X(k)|’);

شکل ۶- پاسخ اندازه نتیجه FFT بعد از اعمال fftshift، رسم شده بر حسب اندیس (بالا) و بر حسب فرکانس محاسبه شده (پایین).

با مقایسه شکل پایینی در دو شکل ۴ و ۶ میتوان دید که محور فرکانسی مرتب شده برای تفسیر محتوای فرکانسی مناسب‌تر است.

IFFTshift

اثر تابع fftshift با استفاده تابع معکوس آن، ifftshift قابل برگشت است. تابع ifftshift ترتیب فرکانسی خام اولیه را بازیابی میکند. اگر خروجی FFT بوسیله تابع fftshift مرتب شده باشد، آنگاه حتما قبل از اعمال تبدیل معکوس فوریه ifft باید اجزاء‌ فرکانسی به ترتیب اولیه خود برگردند.

X = fft(x,N) %compute X[k]

x = ifft(X,N) %compute x[n]

X = fftshift(fft(x,N)); %take FFT and rearrange frequency order (this is mainly done for interpretation)

x = ifft(ifftshift(X),N) % restore raw frequency order and then take IFFT

برخی ملاحظات در رابطه با FFTshift‌ و IFFTshift

زمانیکه N‌ یک عدد فرد است، برای یک دنباله دلخواه، توابع fftshift و ifftshift نتایج متفاوتی تولید میکنند. هر چند که، وقتی این دو تابع با هم استفاده شوند، دنباله اصلی قابل بازیابی است.

>>  x=[0,1,2,3,4,5,6,7,8]
0     1     2     3     4     5     6     7     8
>>  fftshift(x)
5     6     7     8     0     1     2     3     4
>>  ifftshift(x)
4     5     6     7     8     0     1     2     3
>>  ifftshift(fftshift(x))
0     1     2     3     4     5     6     7     8
>>  fftshift(ifftshift(x))
0     1     2     3     4     5     6     7     8

وقتی N‌ عدد زوج است، برای یک دنباله دلخواه، توابع fftshift‌ و ifftshift نتایج یکسانی تولید میکنند. وقتی آنها با هم استفاده شوند، دنباله اصلی قابل بازیابی است.

>>  x=[0,1,2,3,4,5,6,7]
0     1     2     3     4     5     6     7
>>  fftshift(x)
4     5     6     7     0     1     2     3
>>  ifftshift(x)
4     5     6     7     0     1     2     3
>>  ifftshift(fftshift(x))
0     1     2     3     4     5     6     7
>>  fftshift(ifftshift(x))
0     1     2     3     4     5     6     7

منبع: https://www.gaussianwaves.com

تفسیر تبدیل فوریه سریع (FFT)، تبدیل فوریه گسسته (DFT) مختلط، بازه‌های فرکانسی و fftshift

دیدگاه ها (0)

لغو پاسخ