تبدیل فوریه سریع (fft) — به زبان ساده
تبدیل فوریه سریع (Fast Fourier Transform) یا FFT یکی از مهمترین الگوریتمهای مورد استفاده در پردازش سیگنال و آنالیز داده است. در واقع FFT یک الگوریتم است که برای محاسبه تبدیل فوریه گسسته (Discrete Fourier Transform) یا DFT و نیز معکوس آن (IDFT) مورد استفاده قرار میگیرد. در این مطلب قصد داریم به بیان تبدیل فوریه سریع بپردازیم و روشهای محاسبه آن را بررسی کنیم.
آنالیز فوریه میتواند یک سیگنال از حوزه اصلی، که معمولا زمان یا فضا است را به نمایشی در حوزه فرکانس و نیز بلعکس تبدیل کند. تبدیل فوریه گسسته معمولا از طریق تجزیه دنباله مقادیر، به عناصر با فرکانسهای متفاوت محاسبه میشود. این تبدیل در بسیاری از رشتهها مفید است، اما مشکلی که وجود دارد این است که محاسبه مستقیم این تبدیل با استفاده از تعریف آن بسیار کند است و در عمل کاربردی ندارد. تبدیل فوریه سریع یا FFT روشی است که به وسیله آن میتوان تبدیل فوریه گسسته را به سرعت محاسبه کرد. در واقع تبدیل فوریه سریع از طریق تجزیه ماتریس DFT به حاصلضرب ماتریسهای تنک (Sparse) که در آنها اکثر داریههای ماتریس صفر هستند، محاسبات را تسریع میبخشد.
در نتیجه، تبدیل فوریه سریع قادر است که پیچیدگی (Complexity) محاسبات DFT را از $$ \mathcal{O}[N^2] $$ به $$ \mathcal{O}[N\log N] $$ کاهش دهد. تبدیل فوریه گسسته مانند نوع آشناتر خود، یعنی تبدیل فوریه پیوسته، دارای فرم مستقیم و معکوس است.
کاربرد تبدیل فوریه
تبدیل فوریه سریع یا Fast Fourier Transform که به آن fft نیز میگویند، در علوم مهندسی میتواند برای تعیین فرکانسهای غالب یک سیگنال ارتعاشی (Vibration Signal) به کار رود. زمانی که فرکانسهای غالب (Dominant Frequencies) یک سیستم متناظر با فرکانسهای طبیعی (Natural Frequency) آن باشند، ارتعاشهای رخ داده میتوانند به دلیل رزونانس (Resonance) تقویت شوند.
این اتفاق میتواند تا حدی ادامه یابد که منجر به آسیب دیدن یا ریزش و تخریب یک سازه شود. فرض کنید یک سیستم صوتی در اختیار داریم و فرکانس طبیعی پنجره اتاق در حدود ۱۰۰ هرتز است. حال با استفاده از تبدیل فوریه میتوانیم بدانیم که تا چه حد مجاز هستیم صدای صوت را بالا ببریم. در ادامه به بیان روشی برای پاسخ به این سوال و نیز سوالات مشابه میپردازیم.
سیگنال در حوزه زمان
تبدیل فوریه معمولا برای تبدیل یک سیگنال در طیف زمانی به سیگنالی در طیف فرکانسی مورد استفاده قرار میگیرد. مثالهایی از امواج در طیف زمانی را میتوان سیگنالهای صوتی، سیگنالهای الکتریکی و ارتعاشات مکانیکی نام برد. تصویر زیر نمونهای از یک سیگنال حوزه زمان را در طول ۰٫۲۵ ثانیه نشان میدهد.
همانطور که از تصویر دیده میشود، این سیگنال شبیه به یک موج با فرکانسهای مختلف و یا ترکیب چند موج است.
تبدیل فوریه گسسته
حال میخواهیم از تبدیل فوریه استفاده کنیم. تبدیل فوریه در واقع از این اصل استفاده میکند که هر تابع غیر خطی را میتوان توسط مجموع بینهایتِ امواج سینوسی نمایش داد. در تصویر زیر نمایی از مفهوم نمایش یک سیگنال توسط مجموع سیگنالهای سینوسی نمایش داده شده است.
در واقع تصویر بالا نشان میدهد که یک تابع پله چگونه با استفاده از مجموع چندین موج سینوسی بازسازی میشود. نمودار قرمز رنگ نشان دهنده تابع پله در حوزه زمان و نمودار آبی رنگ نشان دهنده این تابع در حوزه فرکانس است. مجموع امواج سنوسی کوچک در پس زمینه تصویر میتوانند تابع پله را بازسازی کنند.
تبدیل فوریه، یک سیگنال زمانی را تجزیه میکند و اطلاعاتی راجع به فرکانسهای تمام امواج سینوسی مورد نیاز برای ساخت سیگنال حوزه زمان را در اختیار ما قرار میدهد.
برای یک دنباله با مقادیر فاصله برابر، فرمول تبدیل فوریه گسسته مستقیم به صورت زیر بیان میشود:
$$ X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~k~n~/~N} $$
همچنین فرمول تبدیل فوریه گسسته معکوس به صورت زیر است:
$$ x_n = \frac{1}{N}\sum_{k=0}^{N-1} X_k e^{i~2\pi~k~n~/~N} $$
در فرمول بالا، $$ N $$ تعداد نمونهها (Sample) و $$ n $$ نمونه فعلی است. همچنین $$ x_n $$ مقدار سیگنال در زمان n و $$ k $$ فرکانس فعلی (۰ هرتز تا N-1 هرتز) و $$ X_k $$ نتیجه حاصل از تبدیل فوریه گسسته است.
در واقع تبدیل از $$ x_n \to X_k $$، یک تبدیل از فضای پیکربندی به فضای فرکانس است. این تبدیل هم در بررسی طیف توان سیگنال و هم در انتقال برخی مسائل خاص به فضایی جهت انجام راحتتر محاسبات، بسیار مفید است.
برای سادگی ما فقط به تبدیل فوریه گسسته مستقیم میپردازیم؛ زیرا تبدیل فوریه گسسته معکوس را نیز میتوان به صورت مشابه محاسبه کرد. با نگاهی دقیقتر به فرمول تبدیل فوریه گسسته در بالا، میتوان به این نکته پی برد که این فرمول چیزی بیش از ضرب نقطهای نیست.
همان طور که به یاد داریم، ضرب نقطهای به صورت زیر ایجاد میشود:
$$ a \cdot b = \sum_{n=1}^{n}a_ib_i $$
بنابراین، فرمول بالا را میتوان به این صورت بازنویسی کرد:
$$ \overrightarrow{X} = M \overrightarrow{x} $$
که در این فرمول ماتریس $$ M $$ به صورت زیر نوشته میشود:
$$ M_{kn} = e^{-i~2\pi~k~n~/~N}. $$
برای محاسبه ضرب نقطهای در پایتون از کد زیر استفاده میشود:
1import numpy as np
2
3def DFT(x):
4 """
5 Compute the discrete Fourier Transform of the 1D array x
6 :param x: (array)
7 """
8
9 N = x.size
10 n = np.arange(N)
11 k = n.reshape((N, 1))
12 e = np.exp(-2j * np.pi * k * n / N)
13 return np.dot(e, x)
اما باید به این نکته توجه کرد که اگر این برنامه را بر روی سیگنال حوزه زمان اجرا کنیم، که حدودا شامل ۱۰۰۰۰ نمونه است، آنگاه اجرای آن حدود ۱۰ ثانیه به طول میانجامد و با توجه به اینکه معمولا این دستور فقط بخشی از برنامه کلی مورد نیاز برای پروژههای مختلف است، میتوان گفت سرعت بسیار پایینی دارد.
پیچیدگی محاسباتی تبدیل فوریه گسسته
همان طور که میدانیم، معادله تبدیل فوریه گسسته برای N نقطه از دنباله با طول محدود $$ x(n) $$ توسط فرمول زیر محاسبه میشود:
$$ X(k)=\sum\limits_{n=0}^{N-1}{x(n){{e}^{-j\tfrac{2\pi }{N}kn}}} $$
حال میخواهیم تعداد ضرب و جمعهای مورد نیاز برای محاسبه تبدیل فوریه گسسته با استفاده از فرمول بالا را به دست آوریم. برای هر ضریب در تبدیل فوریه گسسته $$ X(k) $$، باید N عبارت را که شامل $$ x(1){{e}^{-j\tfrac{2\pi }{N}k \times 1}} $$، $$ x(0){{e}^{-j\tfrac{2\pi }{N}k \times 0}} $$ تا $$ x(N-1){{e}^{-j\tfrac{2\pi }{N}k \times (N-1)}} $$ هستند، محاسبه کنیم. سپس باید این عبارتها را با هم جمع کنیم.
در حالت کلی $$ x(n) $$ یک عدد مختلط است و هر خروجی تبدیل فوریه گسسته تقریبا به N ضرب مختلط و N-1 جمع مختلط نیاز دارد. یک ضرب مختلط به تنهایی به ۴ ضرب حقیقی و ۲ جمع حقیقی نیاز دارد. به عنوان مثال، ضرب عدد مختلط $$ d_1=a_1+jb_1 $$ در $$ d_2=a_2+jb_2 $$، پاسخی برابر با $$ d_1d_2=(a_1a_2-b_1b_2)+j(b_1a_2+a_1b_2) $$ دارد.
علاوه بر این، یک جمع مختلط خود به ۲ جمع حقیقی نیاز دارد، بنابراین هر ضریب در تبدیل فوریه گسسته به $$ 4N $$ ضرب حقیقی و $$ 2N+2(N-1)=4N-2 $$ جمع حقیقی نیاز دارد. برای محاسبه تبدیل فوریه گسسته نهایی، لازم است تا معادله DFT را برای تمام N مقدار از k محاسبه کنیم. بنابراین آنالیز تبدیل فوریه N نقطه، به $$ 4N^2 $$ جمع حقیقی و $$ N(4N-2) $$ جمع حقیقی نیاز دارد. بنابراین، همان طور که قبلا نیز به آن اشاره کردیم، محاسبات مورد نیاز برای یک تبدیل فوریه گسسته N نقطهای، از درجه $$ N^2 $$ خواهد بود. بنابراین با افزایش N، محاسبات مورد نیاز به سرعت افزایش مییابد. در شکل زیر، تعداد محاسبات ضرب حقیقی مورد نیاز بر حسب طول N تبدیل فوریه گسسته ترسیم شده است.
الگوریتمهای زیادی برای بهبود مشکل تعداد بالای محاسبات و سرعت پایین در تبدیل فوریه گسسته به وجود آمدهاند. در ادامه به بررسی یکی از این الگوریتمها میپردازیم.
تبدیل فوریه سریع
افرادی به نامهای کولی و توکی (Cooley and Tukey) توانستند الگوریتمی برای محاسبه تبدیل فوریه سریع یا Fast Fourier Transform به دست بیاورند. در این الگوریتم که مهمترین الگوریتم تبدیل فوریه سریع است، به صورت بازگشتی (Recursively) تبدیل فوریه گسسته را به مسايل کوچکتر میشکند و زمان مورد نیاز برای انجام محاسبات را به مقدار قابل توجهی کاهش میدهد.
همان طور که قبلا هم اشاره شد، پیچیدگی یک مسئله تبدیل فوریه گسسته از درجه $$ \mathcal{O}[N^2] $$ به $$ \mathcal{O}[N\log N] $$ در تبدیل فوریه سریع کاهش مییابد. به این معنی که اگر $$ N= 10^6 $$ المان داشته باشد، انتظار میرود که یک الگوریتم FFT بتواند مسئله را در زمانی برابر با 50 ثانیه حل کند، اما این زمان برای تبدیل فوریه گسسته برابر با 20 ساعت خواهد بود. حال سوالی که پیش میآید این است که تبدیل فوریه سریع چگونه به این سرعت بالا در محاسبات دست مییابد؟
تقارن در تبدیل فوریه گسسته
یکی از مهمترین ابزارها در ساخت یک الگوریتم در برنامهنویسی این است که تقارن را در حل مسئله اعمال کنیم. اگر به صورت تحلیلی بتوان نشان داد که یک قسمت از مسئله به قسمت دیگری از آن مشابه است، میتوان به سادگی آن قسمت را یک بار محاسبه کرد و هزینه محاسباتی را کاهش داد. کولی و توکی نیز از همین روش استفاده کردند و حل مسئله تبدیل فوریه گسسته را سریعتر کردند.
حال میتوان بیان روش را با این پرسش شروع کرد که مقدار $$ X_{N+k} $$ با توجه به فرمول تبدیل فوریه گسسته بالا چقدر است. پس از جایگذاری مقدار به صورت زیر به دست میآید:
$$ \begin{align*}
X_{N + k} &= \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~(N + k)~n~/~N}\\
&= \sum_{n=0}^{N-1} x_n \cdot e^{- i~2\pi~n} \cdot e^{-i~2\pi~k~n~/~N}\\
&= \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~k~n~/~N}
\end{align*} $$
در این فرمول از خاصیت همانندی (Identity) استفاده شده است که برای هر عدد صحیح n صحیح است.
خط آخر در فرمول بالا، مشخصه تقارن بسیار مهمی را در تبدیل فوریه گسسته نشان میدهد که در نتیجه آن داریم:
$$ X_{N+k} = X_k. $$
با بسط این عبارت به سادگی میتوان نوشت که:
$$ X_{k + i \cdot N} = X_k $$
برای هر i صحیح، همان طور که در بالا نشان داده شد، این تقارن صادق است. بنابراین، تقارن را میتوان برای محاسبه تبدیل فوریه گسسته با سرعت بالاتر نیز اعمال کرد.
اعمال تقارن و تبدیل DFT به FFT
کولی و توکی نشان دادند که این امکان وجود دارد که محاسبات لازم برای تبدیل فوریه گسسته را به دو قسمت کوچکتر بشکنیم. با استفاده از تعریف تبدیل فوریه گسسته داریم:
$$ \begin{align}
X_k &= \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~k~n~/~N} \\
&= \sum_{m=0}^{N/2 - 1} x_{2m} \cdot e^{-i~2\pi~k~(2m)~/~N} + \sum_{m=0}^{N/2 - 1} x_{2m + 1} \cdot e^{-i~2\pi~k~(2m + 1)~/~N} \\
&= \sum_{m=0}^{N/2 - 1} x_{2m} \cdot e^{-i~2\pi~k~m~/~(N/2)} + e^{-i~2\pi~k~/~N} \sum_{m=0}^{N/2 - 1} x_{2m + 1} \cdot e^{-i~2\pi~k~m~/~(N/2)}
\end{align} $$
در این حالت، تبدیل فوریه گسسته تکی را به دو عبارت تقسیم کردیم که هر کدام شباهت بسیار زیادی به عبارت تبدیل فوریه اصلی دارند و یکی بر روی اعداد فرد و دیگری بر روی اعداد زوج عمل میکنند. اما تا این قسمت هنوز هیچ توان محاسباتی را کاهش ندادهایم و هر عبارت از $$ (N/2)*N $$ محاسبه تشکیل شده است که در مجموع $$ N^2 $$ محاسبه را شامل میشود.
آنچه موجب کاهش هزینه محاسباتی میشود، استفاده از خاصیت تقارن در هر یک از عبارتهای محاسبه شده است. زیرا بازه $$ k $$ برابر با $$ 0 \le k < N $$ است، در حالی که بازه $$ n $$ برابر با $$ 0 \le n < M \equiv N/2 $$ در نظر گرفته میشود. در خاصیت تقارن بالا نکتهای که وجود دارد این است که میتوان برای هر زیر مسئله (Sub-Problem) فقط نیمی از محاسبات را انجام داد. بنابراین محاسبات $$ \mathcal{O}[N^2] $$ تبدیل به $$ \mathcal{O}[M^2] $$ میشود که $$ M $$ برابر با نصف اندازه $$ N $$ است.
اما روند کار در اینجا متوقف نمیشود. تا زمانی که تبدیل فوریه کوچکتر دارای مقدار M زوج باشد، میتوانیم این روش تقسیم و غلبه (Divide-and-Conquer) را به صورت تکراری انجام دهیم و هر بار هزینه محاسباتی را نصف کنیم. این روند را تا جایی ادامه میدهیم که آرایه به دست آمده آنقدر کوچک باشد که استفاده از این استراتژی دیگر تاثیری در بهبود محاسبات نداشته باشد.
مثالی از تبدیل فوریه سریع با $$ N=8 $$
فرض کنید که $$ x(n) $$ دنبالهای متشکل از ۸ مقدار باشد. با استفاده از فرمول تبدیل فوریه گسسته، میتوانیم مقدار DFT این سیگنال را به صورت زیر به دست آوریم:
$$ X\left(k\right)=x(0){{e}^{-j\tfrac{2\pi }{8}k\times0}}+x(1){{e}^{-j\tfrac{2\pi }{8}k\times1}}+x(2){{e}^{-j\tfrac{2\pi }{8}k\times2}}+\dots +x(7){{e}^{-j\tfrac{2\pi }{8}k\times7}} $$
همان طور که در بالا بیان شد، هدف این است که بتوانیم این دنباله هشت نقطهای را ابتدا توسط دو عبارت با طول کوتاهتر بازنویسی کنیم. طول DFT در این حالت برابر با $$ N=8 $$ است که در مخرج عبارت نمایی مختلط ظاهر میشود. در ایده تبدیل فوریه سریع، میخواهیم برخی از عبارتهای معادله تبدیل فوریه گسسته را به نحوی گروهبندی کنیم که کسر $$ -\tfrac{j2\pi}{8}kn $$ سادهتر شود.
بنابراین، DFTهایی با طول کوتاهتر را از فرمول اصلی استخراج میکنیم. مثلا در این حالت که N عددی زوج است، تمام عبارتها با اندیس نمونه زوج، یعنی $$ x(4) $$ و $$ x(2) $$ و $$ x(0) $$ و $$ x(6) $$ را انتخاب میکنیم. بنابراین عبارت زیر را به دست میآوریم.
$$ G\left(k\right)=x(0){{e}^{-j\tfrac{2\pi }{8}k\times0}}+x(2){{e}^{-j\tfrac{2\pi }{8}k\times2}}+x(4){{e}^{-j\tfrac{2\pi }{8}k\times4}}+x(6){{e}^{-j\tfrac{2\pi }{8}k\times6}} $$
پس از سادهسازی کسرها در عبارات مختلط، مقدار زیر محاسبه میگردد:
$$ G\left(k\right)=x(0){{e}^{-j\tfrac{2\pi }{4}k\times0}}+x(2){{e}^{-j\tfrac{2\pi }{4}k\times1}}+x(4){{e}^{-j\tfrac{2\pi }{4}k\times2}}+x(6){{e}^{-j\tfrac{2\pi }{4}k\times3}} $$
با مقایسه این عبارت با فرمول اصلی محاسبه DFT مشاهده میکنیم که $$ G\left(k\right) $$ یک تبدیل فوریه گسسته ۴ نقطهای متشکل از $$ x(4) $$ و $$ x(2) $$ و $$ x(0) $$ و $$ x(6) $$ است. تا این قسمت، نشان دادیم که نصف محاسبات در فرمول اولیه را میتوانیم با محاسبه ۴ نقطه از تبدیل فوریه گسسته حذف کنیم. اما برای مابقی محاسبات، که متعلق به نمونههای با اندیس فرد هستند، عبارت زیر را داریم:
$$ H_{1}\left(k\right)=x(1){{e}^{-j\tfrac{2\pi }{8}k\times1}}+x(3){{e}^{-j\tfrac{2\pi }{8}k\times3}}+x(5){{e}^{-j\tfrac{2\pi }{8}k\times5}}+x(7){{e}^{-j\tfrac{2\pi }{8}k\times7}} $$
برای سادهسازی $$ -\tfrac{j2\pi}{8}kn $$، میتوانیم از $$ e^{-j\tfrac{2\pi}{8}k} $$ فاکتور بگیریم و عبارت بالا را بازنویسی کنیم:
$$ H_{1}\left(k\right)= {{e}^{-j\tfrac{2\pi }{8}k\times1}}\left( x(1){{e}^{-j\tfrac{2\pi }{8}k\times0}}+x(3){{e}^{-j\tfrac{2\pi }{8}k\times2}}+x(5){{e}^{-j\tfrac{2\pi }{8}k\times4}}+x(7){{e}^{-j\tfrac{2\pi }{8}k\times6}} \right) $$
مقدار بالا را مجددا سادهسازی میکنیم:
$$ H_{1}\left(k\right)= {{e}^{-j\tfrac{2\pi }{8}k\times1}}\left( x(1){{e}^{-j\tfrac{2\pi }{4}k\times0}}+x(3){{e}^{-j\tfrac{2\pi }{4}k\times1}}+x(5){{e}^{-j\tfrac{2\pi }{4}k\times2}}+x(7){{e}^{-j\tfrac{2\pi }{4}k\times3}} \right) $$
حال بار دیگر این فرمول را با فرمول اولیه برای محاسبه DFT مقایسه میکنیم و میبینیم که عبارت از ضرب DFT چهار نقطه یعنی $$ x(۵) $$ و $$ x(۳) $$ و $$ x(۱) $$ و $$ x(۷) $$ به دست آمده است. بنابراین، به هدف محاسبه DFT هشت نقطهای با استفاده از دو DFT کوتاهتر ۴ نقطهای رسیدیم. نتیجه نهایی را میتوان به صورت زیر نمایش داد:
$$ X\left(k\right)=G\left(k\right)+{{e}^{-j\tfrac{2\pi }{8}k\times1}}H\left(k\right) $$
که در این فرمول، $$ G(k) $$ و $$ H(k) $$ متناظر با تبدیل فوریه گسسته ۴ نقطهای هستند. اما سوالی که به وجود میآید این است که چگونه $$ G(k)+e^{-\tfrac{j2\pi}{8}k \times 1}H(k) $$ بار محاسباتی را کاهش میدهد، در حالی که به تعداد برابری ضرب حقیقی در مقایسه با $$ X(k) $$ در فرمول اصلی نیاز دارد؟ در واقع ما فقط فرمول را بازآرایی کردیم و فاکتور ضرایب را تغییر دادیم، اما تعداد ضربها اصلا کاهش نیافتهاند.
برای حل این تضاد ظاهری، باید سه نکته را به یاد داشته باشیم. اولا، یک عبارت نمایی مختلط گسسته $$ e^{-\tfrac{j2\pi}{N}kn} $$ یک تابع متناوب از k با تناوب N است. دوما، بازه k در این مثال از ۰ تا ۷ است.
سوما، زمانی که معادله $$ X\left(k\right)=x(0){{e}^{-j\tfrac{2\pi }{8}k\times0}}+x(1){{e}^{-j\tfrac{2\pi }{8}k\times1}}+x(2){{e}^{-j\tfrac{2\pi }{8}k\times2}}+\dots +x(7){{e}^{-j\tfrac{2\pi }{8}k\times7}} $$ متشکل از عبارات نمایی با تناوب ۸، یعنی $$ e^{-j\tfrac{2\pi}{8}kn} $$ باشد، آنگاه معادله $$ X\left(k\right)=G\left(k\right)+{{e}^{-j\tfrac{2\pi }{8}k\times1}}H\left(k\right) $$ متشکل از دو تبدیل فوریه گسسته ۴ نقطه|ی $$ G(k) $$ و $$ H(k) $$ است که مجموع عبارات نمایی مختلط با تناوب ۴ هستند.
این رفتار متناوب، به ما اجازه میدهد تا پیچیدگی محاسباتی کلی را کاهش دهیم.
به عنوان مثال، زمانی که $$ X(k) $$ و $$ X(k+4) $$ را محاسبه میکنیم، بر اساس معادله اول، مجبوریم که عبارات معادله را دو بار حساب کنیم. اما با استفاده از معادله دوم که به دست آوردیم، داریم:
$$ X\left(k\right)=G\left(k\right)+{{e}^{-j\tfrac{2\pi }{8}k\times1}}H\left(k\right) $$
$$ X\left(k+4\right)=G\left(k+4\right)+{{e}^{-j\tfrac{2\pi }{8}(k+4)\times1}}H\left(k+4\right) $$
چون $$ G(k) $$ و $$ H(k) $$ متناوب با دوره تناوب ۴ هستند، بنابراین معادله آخر را میتوان به صورت زیر بازنویسی کرد:
$$ X\left(k+4 \right)=G\left(k\right)+{{e}^{-j\tfrac{2\pi }{8}\left (k+4 \right) \times1}}H\left(k\right) $$
با مقایسه این معادله با معادله اصلی تبدیل فوریه گسسته، مشخص میشود که چگونه شکستن DFT هشت نقطهای به دو DFT چهار نقطهای، باعث میشود از محاسباتی تقریبا یکسان برای $$ X(k) $$ و $$ X(k+4) $$ استفاده کنیم. در نتیجه به صورت قابل توجهی تعداد محاسبات از طریق مشخصه تناوب این توابع، کاهش مییابد.
تصویر زیر نمایی از تجزیه تبدیل فوریه گسسته ۸ نقطهای به دو تبدیل فوریه ۴ نقطهای را نشان میدهد. در این نمودار $$ W_{N}^{k}=e^{-j\tfrac{2\pi}{N}k} $$ است.
زمانی که از معادله اصلی تبدیل فوریه گسسته استفاده میکنیم، نیاز است تا از $$ 4N^2=256 $$ ضرب حقیقی برای هشت نقطه استفاده کنیم، در حالی که با الگوریتم توضیح داده شده، میتوان دو تبدیل فوریه ۴ نقطهای را با استفاده از N ضرب مختلط و $$ 2\times 4 \times (\tfrac{N}{2})^2=2\times 4 \times 4^2=128 $$ ضرب حقیقی محاسبه کرد. N ضرب مختلط خود به $$ 4N=4 \times 8=32 $$ ضرب حقیقی نیاز دارند. بنابراین الگوریتم در حالت کلی به ۱۶۰ ضرب حقیقی نیاز دارد. این عدد در مقایسه با 256 ضرب حقیقی در حالت قبلی، کمتر است. این بهبود قابل توجهی است، اما میتوان از طریق اعمال روش یاد شده به DFT چهار نقطهای به بهبود بیشتری نیز دست یافت. دیدیم که نقطه شروع الگوریتم، طول DFT است که زوج در نظر گرفته میشود و ما قادریم که عبارت نمایی مختلط از اندیس زوج را ساده کنیم. به همین دلیل، برای اینکه بار دیگر الگوریتم را سادهسازی کنیم، نیاز داریم که $$ \tfrac{N}{2} $$ مضربی از دو باشد. در این مثال خاص، $$ \tfrac{N}{2}=4 $$ است و با استفاده از الگوریتم بالا میتوانیم هر کدام از دو تبدیل فوریه ۴ نقطهای را به دو تبدیل فوریه ۲ نقطهای تجزیه کنیم. به عنوان مثال، تبدیل فوریه گسسته $$ \tfrac{N}{2} $$ نقطه اولی در نمودار شکل بالا را میتوان به صورت زیر نشان داد.
با جایگذاری نمودار بالا در هر کدام از تبدیل فوریههای گسسته ۴ نقطهای، میتوان به نمودار جدید زیر رسید.
توجه داشته باشید که این نمودار، هر کدام از ضرایب $$ W_{\tfrac{N}{2}} $$ را با ضریب $$ W_{N}^{2} $$ جایگزین میکند؛ زیرا $$ e^{-j\tfrac{2\pi}{(\tfrac{N}{2})}}=e^{-j\tfrac{2\pi}{N} \times 2} $$ است.
سرانجام، میتوانیم تبدیل فوریه دو نقطهای در تصویر بالا را نیز با نمودار زیر جایگزین کنیم.
بعد از جایگذاری، در نهایت به نمودار زیر میرسیم.
در مثال بررسی شده، نشان داده شد که به منظور پیادهسازی تکنیک به صورت تکراری، تا زمان باقی ماندن تبدیل فوریه گسسته دو نقطهای، باید طول تبدیل فوریه گسسته اصلی را از توان دو انتخاب کنیم. به عبارت دیگر باید $$ N=2^{\nu} $$ برقرار باشد. در این حالت پس از $$ log_{2}(N) $$ تکرار، به دو تبدیل فوریه گسسته دو نقطهای میرسیم. همان طور که از تصویر آخر نیز میتوان مشاهده کرد، هر گام از الگوریتم به $$ Nlog_{2}(N) $$ ضرب مختلط نیاز دارد. مثلا در مورد تبدیل فوریه هشت نقطهای $$ 8log_{2}(8)=24 $$ ضرب مختلط و $$ 24 \times 4=96 $$ ضرب حقیقی نیاز داریم، در حالی که محاسبه مستقیم این DFT به ۲۵۶ ضرب حقیقی نیاز دارد. بنابراین در محاسبه تبدیل فوریه گسسته برای دادههای طولانی، این الگوریتم میتواند منجر به کاهش محاسبات به اندازه قابل توجهی شود.
پیادهسازی بازگشتی تبدیل فوریه سریع
این روش بازگشتی را میتوان در پایتون به سادگی پیادهسازی کرد. همان طور که در قسمت قبل هم بیان شد، تبدیل فوریه گسسته را میتوان به صورت زیر در پایتون اجرا کرد:
1import numpy as np
2def DFT_slow(x):
3 """Compute the discrete Fourier Transform of the 1D array x"""
4 x = np.asarray(x, dtype=float)
5 N = x.shape[0]
6 n = np.arange(N)
7 k = n.reshape((N, 1))
8 M = np.exp(-2j * np.pi * k * n / N)
9 return np.dot(M, x)
حال تبدیل فوریه سریع به صورت زیر نوشته میشود:
1def FFT(x):
2 """A recursive implementation of the 1D Cooley-Tukey FFT"""
3 x = np.asarray(x, dtype=float)
4 N = x.shape[0]
5
6 if N % 2 > 0:
7 raise ValueError("size of x must be a power of 2")
8 elif N <= 32: # this cutoff should be optimized
9 return DFT_slow(x)
10 else:
11 X_even = FFT(x[::2])
12 X_odd = FFT(x[1::2])
13 factor = np.exp(-2j * np.pi * np.arange(N) / N)
14 return np.concatenate([X_even + factor[:N / 2] * X_odd,
15 X_even + factor[N / 2:] * X_odd])
با استفاده از قطعه کد زیر میتوان از صحیح بودن برنامه نوشته شده اطمینان حاصل کرد:
1x = np.random.random(1024)
2np.allclose(FFT(x), np.fft.fft(x))
خروجی این کد True است.
میتوانیم زمان اجرای این برنامه را توسط کد زیر محاسبه کرد:
1%timeit DFT_slow(x)
2%timeit FFT(x)
3%timeit np.fft.fft(x)
خروجی کد به صورت زیر است:
10 loops, best of 3: 77.6 ms per loop 100 loops, best of 3: 4.07 ms per loop 10000 loops, best of 3: 24.7 µs per loop
این روش دارای سرعت بسیار بالاتری نسبت به ورژن اصلی تبدیل فوریه گسسته است. الگوریتم بازگشتی به صورت مجانبی $$ \mathcal{O}[N\log N] $$ است. آنچه که پیادهسازی شد، تبدیل فوریه سریع نام دارد.
به این نکته توجه کنید که هنوز با سرعت الگوریتم توکار (Built-in) برای FFT در پکیج numpy فاصله داریم و البته این اختلاف امری قابل پیش بینی است. الگوریتم FFTPACK در پکیج numpy یک پیادهسازی در فورترن (Fortran) است که حاصل سالها تنظیم و بهینهسازی است. علاوه بر این، از آنجایی که برنامه پیشنهادی ما (که مبتنی بر کتابخانه Numpy است.) از فرایندهایی بازگشتی پشتهای (Python-Stack Recursion) (مبتنی بر پشته) و اختصاص حافظه به تعداد زیادی از آرایههای موقت (Temporary Array) در زبان پایتون استفاده میکند، بار محاسباتی قابل توجهی به سیستم تحمیل میشود.
یک استراتژی مناسب برای افزایش سرعت اجرای کد نوشته شده برای تبدیل فوریه سریع، استفاده از محاسبات برداری تا حد امکان است. برای اجرای کد بالا میتوانیم از این روش استفاده کنیم و فراخوانی توابع بازگشتی را حذف کنیم و برنامه FFT در پایتون را حتی از این هم سریعتر و کاراتر کنیم.
پیادهسازی تبدیل فوریه سریع با استفاده از محاسبات برداری
به این نکته توجه کنید که در پیادهسازی FFT به روش بازگشتی، در پایینترین سطح بازگشت، حاصلضرب بردار در ماتریس را انجام دادیم. کارایی الگوریتم در این حالت میتواند با استفاده از محاسبه ضرب بردار در ماتریس به صورت ضرب یکباره ماتریس در ماتریس بهبود یابد. همچنین، در هر مرحله از برنامه، عملیات تکراری انجام میگیرد که میتواند توسط محاسبات برداری انجام شود. پکیج numpy در این نوع محاسبات بسیار عالی عمل میکند و ما نیز میتوانیم از این ویژگی استفاده کنیم و ورژن برداری برنامه تبدیل فوریه سریع را پیادهسازی کنیم.
1def FFT_vectorized(x):
2 """A vectorized, non-recursive version of the Cooley-Tukey FFT"""
3 x = np.asarray(x, dtype=float)
4 N = x.shape[0]
5
6 if np.log2(N) % 1 > 0:
7 raise ValueError("size of x must be a power of 2")
8
9 # N_min here is equivalent to the stopping condition above,
10 # and should be a power of 2
11 N_min = min(N, 32)
12
13 # Perform an O[N^2] DFT on all length-N_min sub-problems at once
14 n = np.arange(N_min)
15 k = n[:, None]
16 M = np.exp(-2j * np.pi * n * k / N_min)
17 X = np.dot(M, x.reshape((N_min, -1)))
18
19 # build-up each level of the recursive calculation all at once
20 while X.shape[0] < N:
21 X_even = X[:, :X.shape[1] / 2]
22 X_odd = X[:, X.shape[1] / 2:]
23 factor = np.exp(-1j * np.pi * np.arange(X.shape[0])
24 / X.shape[0])[:, None]
25 X = np.vstack([X_even + factor * X_odd,
26 X_even - factor * X_odd])
27
28 return X.ravel()
برنامه در این حالت اندکی مبهم است، اما در واقع یک بازآرایی مجدد از عملگرهای استفاده شده در ورژن بازگشتی برنامه است، با این تفاوت که این بار از تقارن در محاسبه factor استفاده کردهایم و تنها نصفی از آرایه را محاسبه کردهایم. مجددا برای اطمینان از صحت کارکرد برنامه کد زیر را اجرا میکنیم:
1x = np.random.random(1024)
2np.allclose(FFT_vectorized(x), np.fft.fft(x))
خروجی این قطعه کد نیز True است. برای اثبات این مسئله که الگوریتم دارای کارکرد بهینهتری است، میتوان با استفاده از کد زیر بر روی یک آرایه بسیار بزرگ، زمان اجرای این برنامه با برنامههای حالتهای قبل را مقایسه کرد:
1x = np.random.random(1024 * 16)
2%timeit FFT(x)
3%timeit FFT_vectorized(x)
4%timeit np.fft.fft(x)
خروجی این کد به صورت زیر خواهد بود:
10 loops, best of 3: 72.8 ms per loop 100 loops, best of 3: 4.11 ms per loop 1000 loops, best of 3: 505 µs per loop
همان طور که دیده میشود، خروجی این برنامه نسبت به برنامه قبل بهبود یافته است، اما باز هم نسبت به تابع توکار numpy زمان زیادی را برای محاسبات لازم دارد.
طیف فرکانسی تبدیل فوریه سریع
برای درک تاثیر استفاده از تبدیل فوریه سریع، ابتدا کد زیر را برای ایجاد یک سیگنال فرضی در پایتون اجرا میکنیم. توسط این کد ابتدا سیگنالی را تعریف میکنیم که حاصل جمع دو سیگنال سینوسی باشد. یکی از این سیگنالها فرکانس ۴۰ هرتز و دیگری فرکانس ۹۰ هرتز دارد.
1t = np.linspace(0, 0.5, 500)
2s = np.sin(40 * 2 * np.pi * t) + 0.5 * np.sin(90 * 2 * np.pi * t)
3
4plt.ylabel("Amplitude")
5plt.xlabel("Time [s]")
6plt.plot(t, s)
7plt.show()
تصویر زیر، سیگنال را در حوزه زمان نشان میدهد.
برای به دست آوردن طیف فرکانسی سیگنال زمانی به دست آمده در بالا، باید از سیگنال تبدیل فوریه سریع بگیریم. قطعه کد زیر این کار را با استفاده از تابع توکار پایتون انجام میدهد.
1fft = np.fft.fft(s)
2
3
4for i in range(2):
5 print("Value at index {}:\t{}".format(i, fft[i + 1]), "\nValue at index {}:\t{}".format(fft.size -1 - i, fft[-1 - i]))
خروجی این کد به صورت زیر خواهد بود.
Value at index 0: (0.0003804834928402556-0.060555031761900024j) Value at index 499: (0.0003804834928403944+0.060555031761903175j) Value at index 1: (0.0015317714831371565-0.12188808528069561j) Value at index 498: (0.0015317714831373785+0.1218880852806919j)
در قطعه کد بالا، تبدیل فوریه سریع دو سیگنال سینوسی با فرکانس های متفاوت را به دست آوردیم و سپس دو مقدار اول و دو مقدار آخر دنباله تبدیل فوریه سریع را در خروجی نشان دادیم. همان طور که دیده میشود، این مقادیر اعداد مختلط هستند. اگر مقدار اول در این دنباله، که با اندیس ۰ نشان داده شده است را با آخرین عدد از دنباله، که با اندیس ۴۹۹ نشان داده شده است، مقایسه کنیم، میتوانیم به این نکته پی ببریم که قسمت حقیقی هر دو عدد با یکدیگر برابر هستند و مقدار موهومی این اعداد هم دامنه برابری دارند و فقط در علامت متفاوت هستند، یکی از آنها مثبت و دیگری منفی است. در واقع اعداد مختلط مزدوج هستند. این قاعده برای تمام اعداد در دنباله FFT صادق است.
به دلیل این که نیمه دوم دنباله، هیچ اطلاعات جدیدی را به ما نمیدهد، پس میتوانیم نتیجه بگیریم که نیمه اول دنباله FFT، همان خروجی مورد نظر ما محسوب میشود. اعداد خروجی مختلط حاوی اطلاعات زیر هستند:
- دامنه مربوط به فرکانسهای خاصی از موج سینوسی که نشان دهنده انرژی سیگنال هستند.
- انحراف فاز مربوط به فرکانسهای خاصی از موج سینوسی.
دامنه را میتوان با محاسبه قدر مطلق اعداد خروجی و انحراف فاز را با محاسبه زاویه اعداد به دست آورد.
انرژی هر فرکانس از اهمیت بالایی برخوردار است، بنابراین میتوان مقدار مطلق خروجی تبدیل فوریه سریع را محاسبه کرد. برای به دست آوردن درک بهتر از طیف، باید انرژی را بر حسب فرکانس ترسیم کنیم. هر خروجی عدد گسسته FFT، متناظر با یک فرکانس خاص است. رزولوشن فرکانسی توسط فرمول زیر محاسبه میشود:
$$ \Delta f = \frac{f_s}{N} $$
با کنار هم گذاشتن تمام موارد، میتوانیم طیف فرکانسی در مثال ساده سیگنال سینوسی را به دست آوریم. لازم است تا فقط نصفی از طیف را ترسیم کنیم؛ زیرا تمام اطلاعات در همین نیمه به دست میآید. قطعه کد زیر را در پایتون اجرا میکنیم:
1fft = np.fft.fft(s)
2T = t[1] - t[0] # sampling interval
3N = s.size
4
5# 1/T = frequency
6f = np.linspace(0, 1 / T, N)
7
8plt.ylabel("Amplitude")
9plt.xlabel("Frequency [Hz]")
10plt.bar(f[:N // 2], np.abs(fft)[:N // 2] * 1 / N, width=1.5) # 1 / N is a normalization factor
11plt.show()
تصویر زیر به عنوان خروجی این کد ترسیم میشود.
همان طور که میبینیم، تبدیل FFT سیگنال به دست میآید و اطلاعاتی راجع به فرکانسهای امواجی که در حوزه زمان بودند، در اختیار ما قرار میدهد. FFT یک مبادله بین اطلاعات در حوزه زمان و اطلاعات در حوزه فرکانس است. با گرفتن تبدیل FFT از یک سیگنال در حوزه زمان، تمام اطلاعات زمانی سیگنال را از دست داده و در عوض اطلاعات فرکانسی آن را به دست میآوریم. برای داشتن اطلاعات هم در حوزه زمان و هم در حوزه فرکانس در یک طیف، باید از طیف نگاره (Spectrogram) استفاده کنیم.
تصویر زیر FFT یک سیگنال موسیقی را نشان میدهد.
فرکانس این نمودار در مقیاس لگاریتمی ترسیم شده است. همان طور که در تصویر دیده میشود فرکانس غالب این سیگنال بین $$ 10^{1.5} $$ و $$ 10^{2.2} $$ هرتز (30 تا 158 هرتز) است. اگر مطابق مثال ابتدای متن، پنجرهای با فرکانس طبیعی 100 هرتز داشته باشیم، آنگاه میتوان به این نکته پی برد که فرکانس طبیعی پنجره دقیقا در میان این بازه فرکانسی غالب سیگنال صوتی است و با بالا بردن صدای موسیقی، این امکان وجود دارد که پدیده رزونانس به وجود بیاید.
برنامه تبدیل فوریه سریع در جاوا
با استفاده از برنامه زیر میتوان تبدیل فوریه سریع را در زبان برنامهنویسی جاوا محاسبه کرد.
1/*
2 * Free FFT and convolution (Java)
3
4
5
6public final class Fft {
7
8 /*
9 * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
10 * The vector can have any length. This is a wrapper function.
11 */
12 public static void transform(double[] real, double[] imag) {
13 int n = real.length;
14 if (n != imag.length)
15 throw new IllegalArgumentException("Mismatched lengths");
16 if (n == 0)
17 return;
18 else if ((n & (n - 1)) == 0) // Is power of 2
19 transformRadix2(real, imag);
20 else // More complicated algorithm for arbitrary sizes
21 transformBluestein(real, imag);
22 }
23
24
25 /*
26 * Computes the inverse discrete Fourier transform (IDFT) of the given complex vector, storing the result back into the vector.
27 * The vector can have any length. This is a wrapper function. This transform does not perform scaling, so the inverse is not a true inverse.
28 */
29 public static void inverseTransform(double[] real, double[] imag) {
30 transform(imag, real);
31 }
32
33
34 /*
35 * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
36 * The vector's length must be a power of 2. Uses the Cooley-Tukey decimation-in-time radix-2 algorithm.
37 */
38 public static void transformRadix2(double[] real, double[] imag) {
39 // Length variables
40 int n = real.length;
41 if (n != imag.length)
42 throw new IllegalArgumentException("Mismatched lengths");
43 int levels = 31 - Integer.numberOfLeadingZeros(n); // Equal to floor(log2(n))
44 if (1 << levels != n)
45 throw new IllegalArgumentException("Length is not a power of 2");
46
47 // Trigonometric tables
48 double[] cosTable = new double[n / 2];
49 double[] sinTable = new double[n / 2];
50 for (int i = 0; i < n / 2; i++) {
51 cosTable[i] = Math.cos(2 * Math.PI * i / n);
52 sinTable[i] = Math.sin(2 * Math.PI * i / n);
53 }
54
55 // Bit-reversed addressing permutation
56 for (int i = 0; i < n; i++) {
57 int j = Integer.reverse(i) >>> (32 - levels);
58 if (j > i) {
59 double temp = real[i];
60 real[i] = real[j];
61 real[j] = temp;
62 temp = imag[i];
63 imag[i] = imag[j];
64 imag[j] = temp;
65 }
66 }
67
68 // Cooley-Tukey decimation-in-time radix-2 FFT
69 for (int size = 2; size <= n; size *= 2) {
70 int halfsize = size / 2;
71 int tablestep = n / size;
72 for (int i = 0; i < n; i += size) {
73 for (int j = i, k = 0; j < i + halfsize; j++, k += tablestep) {
74 int l = j + halfsize;
75 double tpre = real[l] * cosTable[k] + imag[l] * sinTable[k];
76 double tpim = -real[l] * sinTable[k] + imag[l] * cosTable[k];
77 real[l] = real[j] - tpre;
78 imag[l] = imag[j] - tpim;
79 real[j] += tpre;
80 imag[j] += tpim;
81 }
82 }
83 if (size == n) // Prevent overflow in 'size *= 2'
84 break;
85 }
86 }
87
88
89 /*
90 * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
91 * The vector can have any length. This requires the convolution function, which in turn requires the radix-2 FFT function.
92 * Uses Bluestein's chirp z-transform algorithm.
93 */
94 public static void transformBluestein(double[] real, double[] imag) {
95 // Find a power-of-2 convolution length m such that m >= n * 2 + 1
96 int n = real.length;
97 if (n != imag.length)
98 throw new IllegalArgumentException("Mismatched lengths");
99 if (n >= 0x20000000)
100 throw new IllegalArgumentException("Array too large");
101 int m = Integer.highestOneBit(n) * 4;
102
103 // Trignometric tables
104 double[] cosTable = new double[n];
105 double[] sinTable = new double[n];
106 for (int i = 0; i < n; i++) {
107 int j = (int)((long)i * i % (n * 2)); // This is more accurate than j = i * i
108 cosTable[i] = Math.cos(Math.PI * j / n);
109 sinTable[i] = Math.sin(Math.PI * j / n);
110 }
111
112 // Temporary vectors and preprocessing
113 double[] areal = new double[m];
114 double[] aimag = new double[m];
115 for (int i = 0; i < n; i++) {
116 areal[i] = real[i] * cosTable[i] + imag[i] * sinTable[i];
117 aimag[i] = -real[i] * sinTable[i] + imag[i] * cosTable[i];
118 }
119 double[] breal = new double[m];
120 double[] bimag = new double[m];
121 breal[0] = cosTable[0];
122 bimag[0] = sinTable[0];
123 for (int i = 1; i < n; i++) {
124 breal[i] = breal[m - i] = cosTable[i];
125 bimag[i] = bimag[m - i] = sinTable[i];
126 }
127
128 // Convolution
129 double[] creal = new double[m];
130 double[] cimag = new double[m];
131 convolve(areal, aimag, breal, bimag, creal, cimag);
132
133 // Postprocessing
134 for (int i = 0; i < n; i++) {
135 real[i] = creal[i] * cosTable[i] + cimag[i] * sinTable[i];
136 imag[i] = -creal[i] * sinTable[i] + cimag[i] * cosTable[i];
137 }
138 }
139
140
141 /*
142 * Computes the circular convolution of the given real vectors. Each vector's length must be the same.
143 */
144 public static void convolve(double[] x, double[] y, double[] out) {
145 int n = x.length;
146 if (n != y.length || n != out.length)
147 throw new IllegalArgumentException("Mismatched lengths");
148 convolve(x, new double[n], y, new double[n], out, new double[n]);
149 }
150
151
152 /*
153 * Computes the circular convolution of the given complex vectors. Each vector's length must be the same.
154 */
155 public static void convolve(double[] xreal, double[] ximag,
156 double[] yreal, double[] yimag, double[] outreal, double[] outimag) {
157
158 int n = xreal.length;
159 if (n != ximag.length || n != yreal.length || n != yimag.length
160 || n != outreal.length || n != outimag.length)
161 throw new IllegalArgumentException("Mismatched lengths");
162
163 xreal = xreal.clone();
164 ximag = ximag.clone();
165 yreal = yreal.clone();
166 yimag = yimag.clone();
167 transform(xreal, ximag);
168 transform(yreal, yimag);
169
170 for (int i = 0; i < n; i++) {
171 double temp = xreal[i] * yreal[i] - ximag[i] * yimag[i];
172 ximag[i] = ximag[i] * yreal[i] + xreal[i] * yimag[i];
173 xreal[i] = temp;
174 }
175 inverseTransform(xreal, ximag);
176
177 for (int i = 0; i < n; i++) { // Scaling (because this FFT implementation omits it)
178 outreal[i] = xreal[i] / n;
179 outimag[i] = ximag[i] / n;
180 }
181 }
182
183}
تبدیل فوریه سریع در زبان متلب
با استفاده از برنامه زیر میتوان تبدیل فوریه سریع را در زبان متلب محاسبه کرد.
1function [spectrum, freq] = autofft(xs, ts, fftset)
2
3%% nargin check
4if nargin < 2
5 error('Not enough input arguments.');
6elseif nargin > 3
7 error('Too many input arguments.');
8end
9%
10%% Convert row vectors to column vectors if needed
11if size(xs, 1) == 1 % samples
12 xs = xs(:);
13end
14if size(ts(:), 1) == 1 % sampling frequency
15 fs = ts;
16else
17 fs = 1 / (ts(2) - ts(1));
18end
19%
20%% Specify default fftset
21defset = struct('nwin', size(xs, 1), ...
22 'twin', size(xs, 1) / fs, ...
23 'overlap', 50, ...
24 'lowpass', fs/2, ...
25 'window', 'u', ...
26 'averaging', 'lin', ...
27 'jw', '1', ...
28 'unit', 'pow');
29deffields = fieldnames(defset);
30%
31%% Set analyser parameters
32if nargin == 2 % use default fftset
33 fftset = defset;
34else % use user-defined fftset
35 % Check whether there is user-defined 'nwin' or 'twin' parameter
36 if isfield(fftset, 'nwin')
37 fftset.twin = fftset.nwin / fs;
38 elseif isfield(fftset, 'twin')
39 fftset.nwin = round(fftset.twin * fs);
40 end
41 % Set unspecified parameters to default
42 for i = 1:numel(deffields)
43 if isfield(fftset, deffields{i}) == 0
44 fftset.(deffields{i}) = defset.(deffields{i});
45 end
46 end
47end
48% Generate frequency vector
49freq = (fs * (0:(fftset.nwin/2)) / fftset.nwin)';
50% Set allowed frequencies for the low-pass filtering
51freq = freq(freq <= fftset.lowpass);
52% Set handling of the last spectral line
53if freq(end) == fs/2
54 sh = 1; % Nyquist frequency is at the last spectral line
55else
56 sh = 0; % Nyquist frequency is not at the last spectral line
57end
58% Set number of overlaping samples
59fftset.overlap = round(fftset.nwin * fftset.overlap / 100);
60% Set indices for the signal segmentation
61imax = floor((size(xs, 1)-fftset.overlap) / (fftset.nwin-fftset.overlap));
62 % number of windows
63ind = zeros(imax,2); % matrix of indices
64ni = 1; % pointer
65for i = 1:imax % cycle through windows
66 ind(i,1) = ni;
67 ni = ni + fftset.nwin - 1;
68 ind(i,2) = ni;
69 ni = ni - fftset.overlap + 1;
70end
71% Generate the time weighting window
72[fftset.window, fftset.noiseband] = timeweight(fftset.window, fftset.nwin, fs);
73% Set constant for the jw weigthing
74switch lower(fftset.jw)
75 case '1/jw2'
76 fftset.jw = - 1 ./ (4 * pi^2 * freq.^2);
77 case '1/jw'
78 fftset.jw = 1 ./ (2i * pi * freq);
79 case 'jw'
80 fftset.jw = 2i * pi * freq;
81 case 'jw2'
82 fftset.jw = - 4 * pi^2 * freq.^2;
83 otherwise
84 fftset.jw = 1;
85end
86%
87%% Spectral analyser
88spectrum = zeros(size(freq, 1), size(xs, 2));
89%
90for i = 1:size(xs, 2)
91 % Preallocate an array for temporary spectra
92 spectra = zeros(size(freq, 1), imax);
93 % Frequency analysis of individual segments
94 for j = 1:imax
95 % Fast Fourier transformation of the weighted segment
96 fftSamp = fft(fftset.window .* xs(ind(j,1):ind(j,2), i), ...
97 fftset.nwin) / fftset.nwin;
98 % Application of the jw weigthing and the low-pass filtering
99 spectra(:, j) = fftset.jw .* fftSamp(1:size(freq, 1));
100 % Evaluation of spectral unit
101 switch lower(fftset.unit)
102 case 'rms' % Linear spectrum with rms magnitude
103 spectra(2:end-sh,j) = (2/sqrt(2))*abs(spectra(2:end-sh,j));
104 spectra(end-sh+1:end,j) = abs(spectra(end-sh+1:end,j))/sqrt(2);
105 case 'pk' % Linear spectrum with 0-peak magnitude
106 spectra(2:end-sh,j) = 2 * abs(spectra(2:end-sh,j));
107 case 'pp' % Linear spectrum with peak-peak magnitude
108 spectra(2:end-sh,j) = 4 * abs(spectra(2:end-sh,j));
109 spectra(end-sh+1:end,j) = 2 * abs(spectra(end-sh+1:end,j));
110 case 'psd' % Power spectral density
111 spectra(2:end-sh,j) = sqrt(2) * spectra(2:end-sh,j);
112 spectra(:,j) = (1 / fftset.noiseband) * spectra(:,j) ...
113 .* conj(spectra(:,j));
114 case {'rsd','rmspsd'}% Root mean square of PSD
115 spectra(2:end-sh,j) = sqrt(2) * spectra(2:end-sh,j);
116 spectra(:,j) = sqrt((1 / fftset.noiseband) * ...
117 spectra(:,j) .* conj(spectra(:,j)));
118 otherwise % Autospectrum
119 spectra(2:end-sh,j) = sqrt(2) * spectra(2:end-sh,j);
120 spectra(:,j) = spectra(:,j) .* conj(spectra(:,j));
121 end
122 end
123 % Spectral averaging
124 switch lower(fftset.averaging)
125 case {'energy', 'rms'} % Energy averaging
126 spectrum(:, i) = sqrt(sum(spectra.^2, 2) ./ imax);
127 case {'max', 'peak'} % Maximum peak hold averaging
128 spectrum(:, i) = max(spectra, [], 2);
129 case 'min' % Minimum peak hold averaging
130 spectrum(:, i) = min(spectra, [], 2);
131 case 'none' % No averaging
132 if size(xs, 2) == 1
133 spectrum = spectra;
134 else
135 spectrum{i} = spectra;
136 end
137 otherwise % Linear averaging
138 spectrum(:, i) = mean(spectra, 2);
139 end
140 % Remove imaginary part (residue due to spectra .* conj(spectra))
141 spectrum = real(spectrum);
142end
143% End of main fucntion
144end
145%
146%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147%% Subfunction weightsig generates the time weighting window
148%
149% Input
150% - sym - symbol for the time weighting window
151% - n - length of the time weighting window (samples)
152% - fs - sampling frequency (Hz)
153%
154% Output
155% - window - the time weighting window
156% - noiseband - the noise power bandwidth of the time weighting window
157%
158function [window, noiseband] = timeweight(sym, n, fs)
159 % Generate the specified time weighting window
160 switch sym(1)
161 case 'b' % Blackmann-Harris
162 window = blackmanharris(n);
163 case 'f' % flat-top
164 window = flattopwin(n);
165 case 'h' % Hann
166 window = hann(n);
167 case 'k' % Kaiser-Bessel
168 if length(sym) == 1
169 % window with default beta = 0.5
170 window = kaiser(n, 0.5);
171 else
172 % window with user specified beta
173 window = kaiser(n, str2double(sym(2:end)));
174 end
175 case 'm' % Hamming
176 window = hamming(n);
177 otherwise % uniform
178 window = rectwin(n);
179 end
180 % Adjust magnitude of the generated time weighting window
181 window = window / mean(window);
182 % Calculate the noise power bandwidth in Hz
183 noiseband = enbw(window, fs);
184end
تبدیل فوریه سریع در زبان C
با استفاده از برنامه زیر میتوان تبدیل فوریه سریع را در زبان C محاسبه کرد.
1/*
2 * Free FFT and convolution (C)
3 *
4
5#include <math.h>
6#include <stddef.h>
7#include <stdint.h>
8#include <stdlib.h>
9#include <string.h>
10#include "fft.h"
11
12
13// Private function prototypes
14static size_t reverse_bits(size_t x, int n);
15static void *memdup(const void *src, size_t n);
16
17
18bool Fft_transform(double real[], double imag[], size_t n) {
19 if (n == 0)
20 return true;
21 else if ((n & (n - 1)) == 0) // Is power of 2
22 return Fft_transformRadix2(real, imag, n);
23 else // More complicated algorithm for arbitrary sizes
24 return Fft_transformBluestein(real, imag, n);
25}
26
27
28bool Fft_inverseTransform(double real[], double imag[], size_t n) {
29 return Fft_transform(imag, real, n);
30}
31
32
33bool Fft_transformRadix2(double real[], double imag[], size_t n) {
34 // Length variables
35 bool status = false;
36 int levels = 0; // Compute levels = floor(log2(n))
37 for (size_t temp = n; temp > 1U; temp >>= 1)
38 levels++;
39 if ((size_t)1U << levels != n)
40 return false; // n is not a power of 2
41
42 // Trignometric tables
43 if (SIZE_MAX / sizeof(double) < n / 2)
44 return false;
45 size_t size = (n / 2) * sizeof(double);
46 double *cos_table = malloc(size);
47 double *sin_table = malloc(size);
48 if (cos_table == NULL || sin_table == NULL)
49 goto cleanup;
50 for (size_t i = 0; i < n / 2; i++) {
51 cos_table[i] = cos(2 * M_PI * i / n);
52 sin_table[i] = sin(2 * M_PI * i / n);
53 }
54
55 // Bit-reversed addressing permutation
56 for (size_t i = 0; i < n; i++) {
57 size_t j = reverse_bits(i, levels);
58 if (j > i) {
59 double temp = real[i];
60 real[i] = real[j];
61 real[j] = temp;
62 temp = imag[i];
63 imag[i] = imag[j];
64 imag[j] = temp;
65 }
66 }
67
68 // Cooley-Tukey decimation-in-time radix-2 FFT
69 for (size_t size = 2; size <= n; size *= 2) {
70 size_t halfsize = size / 2;
71 size_t tablestep = n / size;
72 for (size_t i = 0; i < n; i += size) {
73 for (size_t j = i, k = 0; j < i + halfsize; j++, k += tablestep) {
74 size_t l = j + halfsize;
75 double tpre = real[l] * cos_table[k] + imag[l] * sin_table[k];
76 double tpim = -real[l] * sin_table[k] + imag[l] * cos_table[k];
77 real[l] = real[j] - tpre;
78 imag[l] = imag[j] - tpim;
79 real[j] += tpre;
80 imag[j] += tpim;
81 }
82 }
83 if (size == n) // Prevent overflow in 'size *= 2'
84 break;
85 }
86 status = true;
87
88cleanup:
89 free(cos_table);
90 free(sin_table);
91 return status;
92}
93
94
95bool Fft_transformBluestein(double real[], double imag[], size_t n) {
96 bool status = false;
97
98 // Find a power-of-2 convolution length m such that m >= n * 2 + 1
99 size_t m = 1;
100 while (m / 2 <= n) {
101 if (m > SIZE_MAX / 2)
102 return false;
103 m *= 2;
104 }
105
106 // Allocate memory
107 if (SIZE_MAX / sizeof(double) < n || SIZE_MAX / sizeof(double) < m)
108 return false;
109 size_t size_n = n * sizeof(double);
110 size_t size_m = m * sizeof(double);
111 double *cos_table = malloc(size_n);
112 double *sin_table = malloc(size_n);
113 double *areal = calloc(m, sizeof(double));
114 double *aimag = calloc(m, sizeof(double));
115 double *breal = calloc(m, sizeof(double));
116 double *bimag = calloc(m, sizeof(double));
117 double *creal = malloc(size_m);
118 double *cimag = malloc(size_m);
119 if (cos_table == NULL || sin_table == NULL
120 || areal == NULL || aimag == NULL
121 || breal == NULL || bimag == NULL
122 || creal == NULL || cimag == NULL)
123 goto cleanup;
124
125 // Trignometric tables
126 for (size_t i = 0; i < n; i++) {
127 unsigned long long temp = (unsigned long long)i * i;
128 temp %= (unsigned long long)n * 2;
129 double angle = M_PI * temp / n;
130 cos_table[i] = cos(angle);
131 sin_table[i] = sin(angle);
132 }
133
134 // Temporary vectors and preprocessing
135 for (size_t i = 0; i < n; i++) {
136 areal[i] = real[i] * cos_table[i] + imag[i] * sin_table[i];
137 aimag[i] = -real[i] * sin_table[i] + imag[i] * cos_table[i];
138 }
139 breal[0] = cos_table[0];
140 bimag[0] = sin_table[0];
141 for (size_t i = 1; i < n; i++) {
142 breal[i] = breal[m - i] = cos_table[i];
143 bimag[i] = bimag[m - i] = sin_table[i];
144 }
145
146 // Convolution
147 if (!Fft_convolveComplex(areal, aimag, breal, bimag, creal, cimag, m))
148 goto cleanup;
149
150 // Postprocessing
151 for (size_t i = 0; i < n; i++) {
152 real[i] = creal[i] * cos_table[i] + cimag[i] * sin_table[i];
153 imag[i] = -creal[i] * sin_table[i] + cimag[i] * cos_table[i];
154 }
155 status = true;
156
157 // Deallocation
158cleanup:
159 free(cimag);
160 free(creal);
161 free(bimag);
162 free(breal);
163 free(aimag);
164 free(areal);
165 free(sin_table);
166 free(cos_table);
167 return status;
168}
169
170
171bool Fft_convolveReal(const double x[], const double y[], double out[], size_t n) {
172 bool status = false;
173 double *ximag = calloc(n, sizeof(double));
174 double *yimag = calloc(n, sizeof(double));
175 double *zimag = calloc(n, sizeof(double));
176 if (ximag == NULL || yimag == NULL || zimag == NULL)
177 goto cleanup;
178
179 status = Fft_convolveComplex(x, ximag, y, yimag, out, zimag, n);
180cleanup:
181 free(zimag);
182 free(yimag);
183 free(ximag);
184 return status;
185}
186
187
188bool Fft_convolveComplex(
189 const double xreal[], const double ximag[],
190 const double yreal[], const double yimag[],
191 double outreal[], double outimag[], size_t n) {
192
193 bool status = false;
194 if (SIZE_MAX / sizeof(double) < n)
195 return false;
196 size_t size = n * sizeof(double);
197
198 double *xr = memdup(xreal, size);
199 double *xi = memdup(ximag, size);
200 double *yr = memdup(yreal, size);
201 double *yi = memdup(yimag, size);
202 if (xr == NULL || xi == NULL || yr == NULL || yi == NULL)
203 goto cleanup;
204
205 if (!Fft_transform(xr, xi, n))
206 goto cleanup;
207 if (!Fft_transform(yr, yi, n))
208 goto cleanup;
209
210 for (size_t i = 0; i < n; i++) {
211 double temp = xr[i] * yr[i] - xi[i] * yi[i];
212 xi[i] = xi[i] * yr[i] + xr[i] * yi[i];
213 xr[i] = temp;
214 }
215 if (!Fft_inverseTransform(xr, xi, n))
216 goto cleanup;
217
218 for (size_t i = 0; i < n; i++) { // Scaling (because this FFT implementation omits it)
219 outreal[i] = xr[i] / n;
220 outimag[i] = xi[i] / n;
221 }
222 status = true;
223
224cleanup:
225 free(yi);
226 free(yr);
227 free(xi);
228 free(xr);
229 return status;
230}
231
232
233static size_t reverse_bits(size_t x, int n) {
234 size_t result = 0;
235 for (int i = 0; i < n; i++, x >>= 1)
236 result = (result << 1) | (x & 1U);
237 return result;
238}
239
240
241static void *memdup(const void *src, size_t n) {
242 void *dest = malloc(n);
243 if (n > 0 && dest != NULL)
244 memcpy(dest, src, n);
245 return dest;
246}
تبدیل فوریه سریع در زبان #C
با استفاده از کد زیر، میتوان تبدیل فوریه سریع را با زبان #C پیادهسازی کرد.
1/*
2 * Free FFT and convolution (C#)
3 *
4
5
6using System;
7using System.Numerics;
8
9
10public sealed class Fft {
11
12 /*
13 * Computes the discrete Fourier transform (DFT) or inverse transform of the given complex vector, storing the result back into the vector.
14 * The vector can have any length. This is a wrapper function. The inverse transform does not perform scaling, so it is not a true inverse.
15 */
16 public static void Transform(Complex[] vector, bool inverse) {
17 int n = vector.Length;
18 if (n == 0)
19 return;
20 else if ((n & (n - 1)) == 0) // Is power of 2
21 TransformRadix2(vector, inverse);
22 else // More complicated algorithm for arbitrary sizes
23 TransformBluestein(vector, inverse);
24 }
25
26
27 /*
28 * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
29 * The vector's length must be a power of 2. Uses the Cooley-Tukey decimation-in-time radix-2 algorithm.
30 */
31 public static void TransformRadix2(Complex[] vector, bool inverse) {
32 // Length variables
33 int n = vector.Length;
34 int levels = 0; // compute levels = floor(log2(n))
35 for (int temp = n; temp > 1; temp >>= 1)
36 levels++;
37 if (1 << levels != n)
38 throw new ArgumentException("Length is not a power of 2");
39
40 // Trigonometric table
41 Complex[] expTable = new Complex[n / 2];
42 double coef = 2 * Math.PI / n * (inverse ? 1 : -1);
43 for (int i = 0; i < n / 2; i++)
44 expTable[i] = Complex.Exp(new Complex(0, i * coef));
45
46 // Bit-reversed addressing permutation
47 for (int i = 0; i < n; i++) {
48 int j = (int)((uint)ReverseBits(i) >> (32 - levels));
49 if (j > i) {
50 Complex temp = vector[i];
51 vector[i] = vector[j];
52 vector[j] = temp;
53 }
54 }
55
56 // Cooley-Tukey decimation-in-time radix-2 FFT
57 for (int size = 2; size <= n; size *= 2) {
58 int halfsize = size / 2;
59 int tablestep = n / size;
60 for (int i = 0; i < n; i += size) {
61 for (int j = i, k = 0; j < i + halfsize; j++, k += tablestep) {
62 Complex temp = vector[j + halfsize] * expTable[k];
63 vector[j + halfsize] = vector[j] - temp;
64 vector[j] += temp;
65 }
66 }
67 if (size == n) // Prevent overflow in 'size *= 2'
68 break;
69 }
70 }
71
72
73 /*
74 * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
75 * The vector can have any length. This requires the convolution function, which in turn requires the radix-2 FFT function.
76 * Uses Bluestein's chirp z-transform algorithm.
77 */
78 public static void TransformBluestein(Complex[] vector, bool inverse) {
79 // Find a power-of-2 convolution length m such that m >= n * 2 + 1
80 int n = vector.Length;
81 if (n >= 0x20000000)
82 throw new ArgumentException("Array too large");
83 int m = 1;
84 while (m < n * 2 + 1)
85 m *= 2;
86
87 // Trignometric table
88 Complex[] expTable = new Complex[n];
89 double coef = Math.PI / n * (inverse ? 1 : -1);
90 for (int i = 0; i < n; i++) {
91 int j = (int)((long)i * i % (n * 2)); // This is more accurate than j = i * i
92 expTable[i] = Complex.Exp(new Complex(0, j * coef));
93 }
94
95 // Temporary vectors and preprocessing
96 Complex[] avector = new Complex[m];
97 for (int i = 0; i < n; i++)
98 avector[i] = vector[i] * expTable[i];
99 Complex[] bvector = new Complex[m];
100 bvector[0] = expTable[0];
101 for (int i = 1; i < n; i++)
102 bvector[i] = bvector[m - i] = Complex.Conjugate(expTable[i]);
103
104 // Convolution
105 Complex[] cvector = new Complex[m];
106 Convolve(avector, bvector, cvector);
107
108 // Postprocessing
109 for (int i = 0; i < n; i++)
110 vector[i] = cvector[i] * expTable[i];
111 }
112
113
114 /*
115 * Computes the circular convolution of the given complex vectors. Each vector's length must be the same.
116 */
117 public static void Convolve(Complex[] xvector, Complex[] yvector, Complex[] outvector) {
118 int n = xvector.Length;
119 if (n != yvector.Length || n != outvector.Length)
120 throw new ArgumentException("Mismatched lengths");
121 xvector = (Complex[])xvector.Clone();
122 yvector = (Complex[])yvector.Clone();
123 Transform(xvector, false);
124 Transform(yvector, false);
125 for (int i = 0; i < n; i++)
126 xvector[i] *= yvector[i];
127 Transform(xvector, true);
128 for (int i = 0; i < n; i++) // Scaling (because this FFT implementation omits it)
129 outvector[i] = xvector[i] / n;
130 }
131
132
133 private static int ReverseBits(int val) {
134 int result = 0;
135 for (int i = 0; i < 32; i++, val >>= 1)
136 result = (result << 1) | (val & 1);
137 return result;
138 }
139
140}
تبدیل فوریه سریع در زبان جاوا اسکریپت
جهت دسترسی به کدهای پیادهسازی تبدیل فوریه سریع یا fft در جاوا اسکریپت، روی این لینک کلیک کنید و فایل مربوطه حاوی کدهای جاوا اسکریپت را دریافت کنید.
اگر مطلب بالا برای شما مفید بوده و علاقهمند به موضوعات مرتبط با آن هستید، پیشنهاد میکنیم آموزشهای زیر را نیز ببینید:
- مجموعه آموزشهای پردازش سیگنال
- آموزش پردازش سیگنالهای دیجیتال با متلب
- مجموعه آموزشهای مهندسی الکترونیک
- آموزش تجزیه و تحلیل سیگنال ها و سیستم ها
- تبدیل موجک — از صفر تا صد
- سری فوریه مختلط — به زبان ساده
- شبکه عصبی در متلب — از صفر تا صد
^^
عالیه؛ واقعا خلاصه، مفید و تا حد زیادی کامل بود.
اگه امکانش هست بحث بیان FFT با ماتریس تنک که در بالا اشاره کردین ولی بهش نپرداختین رو به همراه بحث محاسبه انتگرال های فوریه با استفاده ازFFT (Computing Fourier Integrals Using the FFT) رو هم اضافه کنین. ممنون
ممنون. خدا خیرتان دهد
با درود
بسیار عالی و کارآمد بود.
سپاس فراوان
سلام…خداقوت بابت مطالب خوبتون.
من بامتلب اشنایی زیادی ندارم و برای پایان نامه باید از تبدیل فوریه استفاده کنم ، منتهی داده هایی از خارج محیط متلب رو وارد متلب کردم و نمیدونم چطور دیتاهای مربوطه رو وارد کدهای تبدیل فوریه که در بالا اورده شده ، کنم…ممنون میشم راهنمایی بفرمایید
عالی بود. ممنون.
ممنون از مطلب جامع و پر بهره تون…دست مریزاد گل کاشتید
با سلام
یک سوال در حوزه فوریه دارم. فرض کنیم ما یک مجموعه از اعداد حقیقی داریم اگر ما از این اعداد تبدیل فوریه بگیریم و در یک ضریب مختلط در همان حوزه فرکانس ضرب کنیم که آن ضریب کسر مختلطی باشدو از حاصل ضرب تبدیل فوریه معکوس بگیریم عددهای حاصل شده ممکن است مختلط باشد؟
با سلام و احترام
ممنون از مطالب خوبتون.
در طی مرور مطالب فوق، یک سوال برام پیش آمد در صورت امکان ممنون میشم راهنمایی بفرمایید.
در کدها بفرض مثال در متلب، ورودی های تابع راتوضیح می دید بایستی چی باشند.
function [spectrum, freq] = autofft(xs, ts, fftset
آیا xs ورودی های تابع در حوزه زمان
ts زمان نمونه برداری
fftset هم تابع مدنظره؟
با تشکر
حمید