تبدیل فوریه سریع (fft) – به زبان ساده

۱۴۶۰۹ بازدید
آخرین به‌روزرسانی: ۲۵ اردیبهشت ۱۴۰۲
زمان مطالعه: ۳۱ دقیقه
دانلود PDF مقاله
تبدیل فوریه سریع (fft) – به زبان سادهتبدیل فوریه سریع (fft) – به زبان ساده

تبدیل فوریه سریع (Fast Fourier Transform) یا FFT یکی از مهم‌ترین الگوریتم‌های مورد استفاده در پردازش سیگنال و آنالیز داده است. در واقع FFT یک الگوریتم است که برای محاسبه تبدیل فوریه گسسته (Discrete Fourier Transform) یا DFT و نیز معکوس آن (IDFT) مورد استفاده قرار می‌گیرد. در این مطلب قصد داریم به بیان تبدیل فوریه سریع بپردازیم و روش‌های محاسبه آن را بررسی کنیم.

997696

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

در نتیجه، تبدیل فوریه سریع قادر است که پیچیدگی (Complexity) محاسبات DFT را از O[N2]\mathcal{O}[N^2] به O[NlogN]\mathcal{O}[N\log N] کاهش دهد. تبدیل فوریه گسسته مانند نوع آشناتر خود، یعنی تبدیل فوریه پیوسته، دارای فرم مستقیم و معکوس است.

کاربرد تبدیل فوریه

تبدیل فوریه سریع یا Fast Fourier Transform که به آن fft نیز می‌گویند، در علوم مهندسی می‌تواند برای تعیین فرکانس‌های غالب یک سیگنال ارتعاشی (Vibration Signal) به کار رود. زمانی که فرکانس‌های غالب (Dominant Frequencies) یک سیستم متناظر با فرکانس‌های طبیعی (Natural Frequency) آن باشند، ارتعاش‌های رخ داده می‌توانند به دلیل رزونانس (Resonance) تقویت شوند.

این اتفاق می‌تواند تا حدی ادامه یابد که منجر به آسیب دیدن یا ریزش و تخریب یک سازه شود. فرض کنید یک سیستم صوتی در اختیار داریم و فرکانس طبیعی پنجره اتاق در حدود ۱۰۰ هرتز است. حال با استفاده از تبدیل فوریه می‌توانیم بدانیم که تا چه حد مجاز هستیم صدای صوت را بالا ببریم. در ادامه به بیان روشی برای پاسخ به این سوال و نیز سوالات مشابه می‌پردازیم.

سیگنال در حوزه زمان

تبدیل فوریه معمولا برای تبدیل یک سیگنال در طیف زمانی به سیگنالی در طیف فرکانسی مورد استفاده قرار می‌گیرد. مثال‌هایی از امواج در طیف زمانی را می‌توان سیگنال‌های صوتی، سیگنال‌های الکتریکی و ارتعاشات مکانیکی نام برد. تصویر زیر نمونه‌ای از یک سیگنال حوزه زمان را در طول ۰٫۲۵ ثانیه نشان می‌دهد.

سیگنال در حوزه زمان
سیگنال در حوزه زمان

همان‌طور که از تصویر دیده می‌شود، این سیگنال شبیه به یک موج با فرکانس‌های مختلف و یا ترکیب چند موج است.

تبدیل فوریه گسسته

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

نمایش یک سیگنال توسط مجموع سیگنال‌های سینوسی
نمایش یک سیگنال توسط مجموع سیگنال‌های سینوسی

در واقع تصویر بالا نشان می‌دهد که یک تابع پله چگونه با استفاده از مجموع چندین موج سینوسی بازسازی می‌شود. نمودار قرمز رنگ نشان‌ دهنده تابع پله در حوزه زمان و نمودار آبی رنگ نشان‌ دهنده این تابع در حوزه فرکانس است. مجموع امواج سنوسی کوچک در پس زمینه تصویر می‌توانند تابع پله را بازسازی کنند.

تبدیل فوریه، یک سیگنال زمانی را تجزیه می‌کند و اطلاعاتی راجع به فرکانس‌های تمام امواج سینوسی مورد نیاز برای ساخت سیگنال حوزه زمان را در اختیار ما قرار می‌دهد.

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

Xk=n=0N1xnei 2π k n / NX_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~k~n~/~N}

همچنین فرمول تبدیل فوریه گسسته معکوس به صورت زیر است:

xn=1Nk=0N1Xkei 2π k n / Nx_n = \frac{1}{N}\sum_{k=0}^{N-1} X_k e^{i~2\pi~k~n~/~N}

در فرمول بالا، NN تعداد نمونه‌ها (Sample) و nn نمونه فعلی است. همچنین xnx_n مقدار سیگنال در زمان n و kk فرکانس فعلی (۰ هرتز تا N-1 هرتز) و XkX_k نتیجه حاصل از تبدیل فوریه گسسته است.

در واقع تبدیل از xnXkx_n \to X_k، یک تبدیل از فضای پیکربندی به فضای فرکانس است. این تبدیل هم در بررسی طیف توان سیگنال و هم در انتقال برخی مسائل خاص به فضایی جهت انجام راحت‌تر محاسبات، بسیار مفید است.

برای سادگی ما فقط به تبدیل فوریه گسسته مستقیم می‌پردازیم؛ زیرا تبدیل فوریه گسسته معکوس را نیز می‌توان به صورت مشابه محاسبه کرد. با نگاهی دقیق‌تر به فرمول تبدیل فوریه گسسته در بالا، می‌توان به این نکته پی برد که این فرمول چیزی بیش از ضرب نقطه‌ای نیست.

همان طور که به یاد داریم، ضرب نقطه‌ای به صورت زیر ایجاد می‌شود:

ab=n=1naibia \cdot b = \sum_{n=1}^{n}a_ib_i

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

X=Mx\overrightarrow{X} = M \overrightarrow{x}

که در این فرمول ماتریس MM به صورت زیر نوشته می‌شود:

Mkn=ei 2π k n / N.M_{kn} = e^{-i~2\pi~k~n~/~N}.

برای محاسبه ضرب نقطه‌ای در پایتون از کد زیر استفاده می‌شود:

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

پیچیدگی محاسباتی تبدیل فوریه گسسته

همان طور که می‌دانیم، معادله تبدیل فوریه گسسته برای N نقطه از دنباله با طول محدود x(n)x(n) توسط فرمول زیر محاسبه می‌شود:

X(k)=n=0N1x(n)ej2πNknX(k)=\sum\limits_{n=0}^{N-1}{x(n){{e}^{-j\tfrac{2\pi }{N}kn}}}

حال می‌خواهیم تعداد ضرب و جمع‌های مورد نیاز برای محاسبه تبدیل فوریه گسسته با استفاده از فرمول بالا را به دست آوریم. برای هر ضریب در تبدیل فوریه گسسته X(k)X(k)، باید N عبارت را که شامل x(1)ej2πNk×1x(1){{e}^{-j\tfrac{2\pi }{N}k \times 1}}، x(0)ej2πNk×0x(0){{e}^{-j\tfrac{2\pi }{N}k \times 0}} تا x(N1)ej2πNk×(N1)x(N-1){{e}^{-j\tfrac{2\pi }{N}k \times (N-1)}} هستند، محاسبه کنیم. سپس باید این عبارت‌ها را با هم جمع کنیم.

در حالت کلی x(n)x(n) یک عدد مختلط است و هر خروجی تبدیل فوریه گسسته تقریبا به N ضرب مختلط و N-1 جمع مختلط نیاز دارد. یک ضرب مختلط به تنهایی به ۴ ضرب حقیقی و ۲ جمع حقیقی نیاز دارد. به عنوان مثال، ضرب عدد مختلط d1=a1+jb1d_1=a_1+jb_1 در d2=a2+jb2d_2=a_2+jb_2، پاسخی برابر با d1d2=(a1a2b1b2)+j(b1a2+a1b2)d_1d_2=(a_1a_2-b_1b_2)+j(b_1a_2+a_1b_2) دارد.

علاوه بر این، یک جمع مختلط خود به ۲ جمع حقیقی نیاز دارد، بنابراین هر ضریب در تبدیل فوریه گسسته به 4N4N ضرب حقیقی و 2N+2(N1)=4N22N+2(N-1)=4N-2 جمع حقیقی نیاز دارد. برای محاسبه تبدیل فوریه گسسته نهایی، لازم است تا معادله DFT را برای تمام N مقدار از k محاسبه کنیم. بنابراین آنالیز تبدیل فوریه N نقطه، به 4N24N^2 جمع حقیقی و N(4N2)N(4N-2) جمع حقیقی نیاز دارد. بنابراین، همان طور که قبلا نیز به آن اشاره کردیم، محاسبات مورد نیاز برای یک تبدیل فوریه گسسته N نقطه‌ای، از درجه N2N^2 خواهد بود. بنابراین با افزایش N، محاسبات مورد نیاز به سرعت افزایش می‌یابد. در شکل زیر، تعداد محاسبات ضرب حقیقی مورد نیاز بر حسب طول N تبدیل فوریه گسسته ترسیم شده است.

نمودار محاسبات مورد نیاز برای DFT با افزایش N
نمودار محاسبات مورد نیاز برای DFT با افزایش N

الگوریتم‌های زیادی برای بهبود مشکل تعداد بالای محاسبات و سرعت پایین در تبدیل فوریه گسسته به وجود آمده‌اند. در ادامه به بررسی یکی از این الگوریتم‌ها می‌پردازیم.

تبدیل فوریه سریع

افرادی به نام‌های کولی و توکی (Cooley and Tukey) توانستند الگوریتمی برای محاسبه تبدیل فوریه سریع یا Fast Fourier Transform به دست بیاورند. در این الگوریتم که مهم‌ترین الگوریتم تبدیل فوریه سریع است، به صورت بازگشتی (Recursively) تبدیل فوریه گسسته را به مسايل کوچک‌تر می‌شکند و زمان مورد نیاز برای انجام محاسبات را به مقدار قابل توجهی کاهش می‌دهد.

همان طور که قبلا هم اشاره شد، پیچیدگی یک مسئله تبدیل فوریه گسسته از درجه O[N2]\mathcal{O}[N^2] به O[NlogN]\mathcal{O}[N\log N] در تبدیل فوریه سریع کاهش می‌یابد. به این معنی که اگر N=106N= 10^6 المان داشته باشد، انتظار می‌رود که یک الگوریتم FFT بتواند مسئله را در زمانی برابر با 50 ثانیه حل کند، اما این زمان برای تبدیل فوریه گسسته برابر با 20 ساعت خواهد بود. حال سوالی که پیش می‌آید این است که تبدیل فوریه سریع چگونه به این سرعت بالا در محاسبات دست می‌یابد؟

تقارن در تبدیل فوریه گسسته

یکی از مهم‌ترین ابزارها در ساخت یک الگوریتم در برنامه‌نویسی این است که تقارن را در حل مسئله اعمال کنیم. اگر به صورت تحلیلی بتوان نشان داد که یک قسمت از مسئله به قسمت دیگری از آن مشابه است، می‌توان به سادگی آن قسمت را یک بار محاسبه کرد و هزینه محاسباتی را کاهش داد. کولی و توکی نیز از همین روش استفاده کردند و حل مسئله تبدیل فوریه گسسته را سریع‌تر کردند.

حال می‌توان بیان روش را با این پرسش شروع کرد که مقدار XN+kX_{N+k} با توجه به فرمول تبدیل فوریه گسسته بالا چقدر است. پس از جایگذاری مقدار به صورت زیر به دست می‌آید:

XN+k=n=0N1xnei 2π (N+k) n / N=n=0N1xnei 2π nei 2π k n / N=n=0N1xnei 2π k n / N\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 صحیح است.

خط آخر در فرمول بالا، مشخصه تقارن بسیار مهمی را در تبدیل فوریه گسسته نشان می‌دهد که در نتیجه آن داریم:

XN+k=Xk.X_{N+k} = X_k.

با بسط این عبارت به سادگی می‌توان نوشت که:

Xk+iN=XkX_{k + i \cdot N} = X_k

برای هر i صحیح، همان طور که در بالا نشان داده شد، این تقارن صادق است. بنابراین، تقارن را می‌توان برای محاسبه تبدیل فوریه گسسته با سرعت بالاتر نیز اعمال کرد.

اعمال تقارن و تبدیل DFT به FFT

کولی و توکی نشان دادند که این امکان وجود دارد که محاسبات لازم برای تبدیل فوریه گسسته را به دو قسمت کوچک‌تر بشکنیم. با استفاده از تعریف تبدیل فوریه گسسته داریم:

Xk=n=0N1xnei 2π k n / N=m=0N/21x2mei 2π k (2m) / N+m=0N/21x2m+1ei 2π k (2m+1) / N=m=0N/21x2mei 2π k m / (N/2)+ei 2π k / Nm=0N/21x2m+1ei 2π k m / (N/2)\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)*N محاسبه تشکیل شده است که در مجموع N2N^2 محاسبه را شامل می‌شود.

آن‌چه موجب کاهش هزینه محاسباتی می‌شود، استفاده از خاصیت تقارن در هر یک از عبارت‌های محاسبه شده است. زیرا بازه kk برابر با 0k<N0 \le k < N است، در حالی که بازه nn برابر با 0n<MN/20 \le n < M \equiv N/2 در نظر گرفته می‌شود. در خاصیت تقارن بالا نکته‌ای که وجود دارد این است که می‌توان برای هر زیر مسئله (Sub-Problem) فقط نیمی از محاسبات را انجام داد. بنابراین محاسبات O[N2]\mathcal{O}[N^2] تبدیل به O[M2]\mathcal{O}[M^2] می‌شود که MM برابر با نصف اندازه NN است.

اما روند کار در اینجا متوقف نمی‌شود. تا زمانی که تبدیل فوریه کوچک‌تر دارای مقدار M زوج باشد، می‌توانیم این روش تقسیم و غلبه (Divide-and-Conquer) را به صورت تکراری انجام دهیم و هر بار هزینه محاسباتی را نصف کنیم. این روند را تا جایی ادامه می‌دهیم که آرایه به دست آمده آنقدر کوچک باشد که استفاده از این استراتژی دیگر تاثیری در بهبود محاسبات نداشته باشد.

مثالی از تبدیل فوریه سریع با N=8N=8

فرض کنید که x(n)x(n) دنباله‌ای متشکل از ۸ مقدار باشد. با استفاده از فرمول تبدیل فوریه گسسته، می‌توانیم مقدار DFT این سیگنال را به صورت زیر به دست آوریم:

X(k)=x(0)ej2π8k×0+x(1)ej2π8k×1+x(2)ej2π8k×2++x(7)ej2π8k×7X\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=8N=8 است که در مخرج عبارت نمایی مختلط ظاهر می‌شود. در ایده تبدیل فوریه سریع، می‌خواهیم برخی از عبارت‌های معادله تبدیل فوریه گسسته را به نحوی گروه‌بندی کنیم که کسر j2π8kn-\tfrac{j2\pi}{8}kn ساده‌تر شود.

بنابراین، DFTهایی با طول کوتاه‌تر را از فرمول اصلی استخراج می‌کنیم. مثلا در این حالت که N عددی زوج است، تمام عبارت‌ها با اندیس نمونه زوج، یعنی x(4)x(4) و x(2)x(2) و x(0)x(0) و x(6)x(6) را انتخاب می‌کنیم. بنابراین عبارت زیر را به دست می‌آوریم.

G(k)=x(0)ej2π8k×0+x(2)ej2π8k×2+x(4)ej2π8k×4+x(6)ej2π8k×6G\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(k)=x(0)ej2π4k×0+x(2)ej2π4k×1+x(4)ej2π4k×2+x(6)ej2π4k×3G\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(k)G\left(k\right) یک تبدیل فوریه گسسته ۴ نقطه‌ای متشکل از x(4)x(4) و x(2)x(2) و x(0)x(0) و x(6)x(6) است. تا این قسمت، نشان دادیم که نصف محاسبات در فرمول اولیه را می‌توانیم با محاسبه ۴ نقطه از تبدیل فوریه گسسته حذف کنیم. اما برای مابقی محاسبات، که متعلق به نمونه‌های با اندیس فرد هستند، عبارت زیر را داریم:

H1(k)=x(1)ej2π8k×1+x(3)ej2π8k×3+x(5)ej2π8k×5+x(7)ej2π8k×7H_{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}}

برای ساده‌سازی j2π8kn-\tfrac{j2\pi}{8}kn، می‌توانیم از ej2π8ke^{-j\tfrac{2\pi}{8}k} فاکتور بگیریم و عبارت بالا را بازنویسی کنیم:

H1(k)=ej2π8k×1(x(1)ej2π8k×0+x(3)ej2π8k×2+x(5)ej2π8k×4+x(7)ej2π8k×6)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)

مقدار بالا را مجددا ساده‌سازی می‌کنیم:

H1(k)=ej2π8k×1(x(1)ej2π4k×0+x(3)ej2π4k×1+x(5)ej2π4k×2+x(7)ej2π4k×3)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(۳) و x(۱)x(۱) و x(۷)x(۷) به دست آمده است. بنابراین، به هدف محاسبه DFT هشت نقطه‌ای با استفاده از دو DFT کوتاه‌تر ۴ نقطه‌ای رسیدیم. نتیجه نهایی را می‌توان به صورت زیر نمایش داد:

X(k)=G(k)+ej2π8k×1H(k)X\left(k\right)=G\left(k\right)+{{e}^{-j\tfrac{2\pi }{8}k\times1}}H\left(k\right)

که در این فرمول، G(k)G(k) و H(k)H(k) متناظر با تبدیل فوریه گسسته ۴ نقطه‌ای هستند. اما سوالی که به وجود می‌آید این است که چگونه G(k)+ej2π8k×1H(k)G(k)+e^{-\tfrac{j2\pi}{8}k \times 1}H(k) بار محاسباتی را کاهش می‌دهد، در حالی که به تعداد برابری ضرب حقیقی در مقایسه با X(k)X(k) در فرمول اصلی نیاز دارد؟ در واقع ما فقط فرمول را بازآرایی کردیم و فاکتور ضرایب را تغییر دادیم، اما تعداد ضرب‌ها اصلا کاهش نیافته‌اند.

برای حل این تضاد ظاهری، باید سه نکته را به یاد داشته باشیم. اولا، یک عبارت نمایی مختلط گسسته ej2πNkne^{-\tfrac{j2\pi}{N}kn} یک تابع متناوب از k با تناوب N است. دوما، بازه k در این مثال از ۰ تا ۷ است.

سوما، زمانی که معادله X(k)=x(0)ej2π8k×0+x(1)ej2π8k×1+x(2)ej2π8k×2++x(7)ej2π8k×7X\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}} متشکل از عبارات نمایی با تناوب ۸، یعنی ej2π8kne^{-j\tfrac{2\pi}{8}kn} باشد، آن‌گاه معادله X(k)=G(k)+ej2π8k×1H(k)X\left(k\right)=G\left(k\right)+{{e}^{-j\tfrac{2\pi }{8}k\times1}}H\left(k\right) متشکل از دو تبدیل فوریه گسسته ۴ نقطه|‌ی G(k)G(k) و H(k)H(k) است که مجموع عبارات نمایی مختلط با تناوب ۴ هستند.

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

به عنوان مثال، زمانی که X(k)X(k) و X(k+4)X(k+4) را محاسبه می‌کنیم، بر اساس معادله اول، مجبوریم که عبارات معادله را دو بار حساب کنیم. اما با استفاده از معادله دوم که به دست آوردیم، داریم:

X(k)=G(k)+ej2π8k×1H(k)X\left(k\right)=G\left(k\right)+{{e}^{-j\tfrac{2\pi }{8}k\times1}}H\left(k\right)

X(k+4)=G(k+4)+ej2π8(k+4)×1H(k+4)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)G(k) و H(k)H(k) متناوب با دوره تناوب ۴ هستند، بنابراین معادله آخر را می‌توان به صورت زیر بازنویسی کرد:

X(k+4)=G(k)+ej2π8(k+4)×1H(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) و X(k+4)X(k+4) استفاده کنیم. در نتیجه به صورت قابل توجهی تعداد محاسبات از طریق مشخصه تناوب این توابع، کاهش می‌یابد.

تصویر زیر نمایی از تجزیه تبدیل فوریه گسسته ۸ نقطه‌ای به دو تبدیل فوریه ۴ نقطه‌ای را نشان می‌دهد. در این نمودار WNk=ej2πNkW_{N}^{k}=e^{-j\tfrac{2\pi}{N}k} است.

تجزیه تبدیل فوریه گسسته ۸ نقطه‌ای به دو تبدیل فوریه ۴ نقطه‌ای
تجزیه تبدیل فوریه گسسته ۸ نقطه‌ای به دو تبدیل فوریه ۴ نقطه‌ای

زمانی که از معادله اصلی تبدیل فوریه گسسته استفاده می‌کنیم، نیاز است تا از 4N2=2564N^2=256 ضرب حقیقی برای هشت نقطه استفاده کنیم، در حالی که با الگوریتم توضیح داده شده، می‌توان دو تبدیل فوریه ۴ نقطه‌ای را با استفاده از N ضرب مختلط و 2×4×(N2)2=2×4×42=1282\times 4 \times (\tfrac{N}{2})^2=2\times 4 \times 4^2=128 ضرب حقیقی محاسبه کرد. N ضرب مختلط خود به 4N=4×8=324N=4 \times 8=32 ضرب حقیقی نیاز دارند. بنابراین الگوریتم در حالت کلی به ۱۶۰ ضرب حقیقی نیاز دارد. این عدد در مقایسه با 256 ضرب حقیقی در حالت قبلی، کمتر است. این بهبود قابل توجهی است، اما می‌توان از طریق اعمال روش یاد شده به DFT چهار نقطه‌ای به بهبود بیشتری نیز دست یافت. دیدیم که نقطه شروع الگوریتم، طول DFT است که زوج در نظر گرفته می‌شود و ما قادریم که عبارت نمایی مختلط از اندیس زوج را ساده کنیم. به همین دلیل، برای اینکه بار دیگر الگوریتم را ساده‌سازی کنیم، نیاز داریم که N2\tfrac{N}{2} مضربی از دو باشد. در این مثال خاص، N2=4\tfrac{N}{2}=4 است و با استفاده از الگوریتم بالا می‌توانیم هر کدام از دو تبدیل فوریه ۴ نقطه‌ای را به دو تبدیل فوریه ۲ نقطه‌ای تجزیه کنیم. به عنوان مثال، تبدیل فوریه گسسته N2\tfrac{N}{2} نقطه اولی در نمودار شکل بالا را می‌توان به صورت زیر نشان داد.

شکستن تبدیل فوریه گسسته <span class=N2\tfrac{N}{2} نقطه در نمودار قبل" width="517" height="276">
شکستن تبدیل فوریه گسسته N2\tfrac{N}{2} نقطه در نمودار قبل

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

محاسبه DFT هشت نقطه‌ای با استفاده از چهار DFT دو نقطه‌ای
محاسبه DFT هشت نقطه‌ای با استفاده از چهار DFT دو نقطه‌ای

توجه داشته باشید که این نمودار، هر کدام از ضرایب WN2W_{\tfrac{N}{2}} را با ضریب WN2W_{N}^{2} جایگزین می‌کند؛ زیرا ej2π(N2)=ej2πN×2e^{-j\tfrac{2\pi}{(\tfrac{N}{2})}}=e^{-j\tfrac{2\pi}{N} \times 2} است.

سرانجام، می‌توانیم تبدیل فوریه دو نقطه‌ای در تصویر بالا را نیز با نمودار زیر جایگزین کنیم.

جایگزین تبدیل فوریه دو نقطه‌ای
جایگزین تبدیل فوریه دو نقطه‌ای

بعد از جایگذاری، در نهایت به نمودار زیر می‌رسیم.

محاسبه تبدیل فوریه تبدیل فوریه گسسته با روش FFT
محاسبه تبدیل فوریه گسسته با روش FFT

در مثال بررسی شده، نشان داده شد که به منظور پیاده‌سازی تکنیک به صورت تکراری، تا زمان باقی ماندن تبدیل فوریه گسسته دو نقطه‌ای، باید طول تبدیل فوریه گسسته اصلی را از توان دو انتخاب کنیم. به عبارت دیگر باید N=2νN=2^{\nu} برقرار باشد. در این حالت پس از log2(N)log_{2}(N) تکرار، به دو تبدیل فوریه گسسته دو نقطه‌ای می‌رسیم. همان‌ طور که از تصویر آخر نیز می‌توان مشاهده کرد، هر گام از الگوریتم به Nlog2(N)Nlog_{2}(N) ضرب مختلط نیاز دارد. مثلا در مورد تبدیل فوریه هشت نقطه‌ای 8log2(8)=248log_{2}(8)=24 ضرب مختلط و 24×4=9624 \times 4=96 ضرب حقیقی نیاز داریم، در حالی که محاسبه مستقیم این DFT به ۲۵۶ ضرب حقیقی نیاز دارد. بنابراین در محاسبه تبدیل فوریه گسسته برای داده‌های طولانی، این الگوریتم می‌تواند منجر به کاهش محاسبات به اندازه قابل توجهی شود.

پیاده‌سازی بازگشتی تبدیل فوریه سریع

این روش بازگشتی را می‌توان در پایتون به سادگی پیاده‌سازی کرد. همان طور که در قسمت قبل هم بیان شد، تبدیل فوریه گسسته را می‌توان به صورت زیر در پایتون اجرا کرد:

حال تبدیل فوریه سریع به صورت زیر نوشته می‌شود:

با استفاده از قطعه کد زیر می‌توان از صحیح بودن برنامه نوشته شده اطمینان حاصل کرد:

خروجی این کد True است.

می‌توانیم زمان اجرای این برنامه را توسط کد زیر محاسبه کرد:

خروجی کد به صورت زیر است:

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

این روش دارای سرعت بسیار بالاتری نسبت به ورژن اصلی تبدیل فوریه گسسته است. الگوریتم بازگشتی به صورت مجانبی O[NlogN]\mathcal{O}[N\log N] است. آن‌چه که پیاده‌سازی شد، تبدیل فوریه سریع نام دارد.

به این نکته توجه کنید که هنوز با سرعت الگوریتم توکار (Built-in) برای FFT در پکیج numpy فاصله داریم و البته این اختلاف امری قابل پیش بینی است. الگوریتم FFTPACK در پکیج numpy یک پیاده‌سازی در فورترن (Fortran) است که حاصل سال‌ها تنظیم و بهینه‌سازی است. علاوه بر این، از آنجایی که برنامه پیشنهادی ما (که مبتنی بر کتابخانه Numpy است.) از فرایندهایی بازگشتی پشته‌ای (Python-Stack Recursion) (مبتنی بر پشته) و اختصاص حافظه به تعداد زیادی از آرایه‌های موقت (Temporary Array) در زبان پایتون استفاده می‌کند، بار محاسباتی قابل توجهی به سیستم تحمیل می‌شود.

یک استراتژی مناسب برای افزایش سرعت اجرای کد نوشته شده برای تبدیل فوریه سریع، استفاده از محاسبات برداری تا حد امکان است. برای اجرای کد بالا می‌توانیم از این روش استفاده کنیم و فراخوانی توابع بازگشتی را حذف کنیم و برنامه FFT در پایتون را حتی از این هم سریع‌تر و کاراتر کنیم.

پیاده‌سازی تبدیل فوریه سریع با استفاده از محاسبات برداری

به این نکته توجه کنید که در پیاده‌سازی FFT به روش بازگشتی، در پایین‌ترین سطح بازگشت، حاصلضرب بردار در ماتریس را انجام دادیم. کارایی الگوریتم در این حالت می‌تواند با استفاده از محاسبه ضرب بردار در ماتریس به صورت ضرب یکباره ماتریس در ماتریس بهبود یابد. همچنین، در هر مرحله از برنامه، عملیات تکراری انجام می‌گیرد که می‌تواند توسط محاسبات برداری انجام شود. پکیج numpy در این نوع محاسبات بسیار عالی عمل می‌کند و ما نیز می‌توانیم از این ویژگی استفاده کنیم و ورژن برداری برنامه تبدیل فوریه سریع را پیاده‌سازی کنیم.

برنامه در این حالت اندکی مبهم است، اما در واقع یک بازآرایی مجدد از عملگرهای استفاده شده در ورژن بازگشتی برنامه است، با این تفاوت که این بار از تقارن در محاسبه factor استفاده کرده‌ایم و تنها نصفی از آرایه را محاسبه کرده‌ایم. مجددا برای اطمینان از صحت کارکرد برنامه کد زیر را اجرا می‌کنیم:

خروجی این قطعه کد نیز True است. برای اثبات این مسئله که الگوریتم دارای کارکرد بهینه‌تری است، می‌توان با استفاده از کد زیر بر روی یک آرایه بسیار بزرگ، زمان اجرای این برنامه با برنامه‌های حالت‌های قبل را مقایسه کرد:

خروجی این کد به صورت زیر خواهد بود:

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 زمان زیادی را برای محاسبات لازم دارد.

طیف فرکانسی تبدیل فوریه سریع

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

تصویر زیر، سیگنال را در حوزه زمان نشان می‌دهد.

سیگنال در حوزه زمان
سیگنال در حوزه زمان

برای به دست آوردن طیف فرکانسی سیگنال زمانی به دست آمده در بالا، باید از سیگنال تبدیل فوریه سریع بگیریم. قطعه کد زیر این کار را با استفاده از تابع توکار پایتون انجام می‌دهد.

خروجی این کد به صورت زیر خواهد بود.

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، متناظر با یک فرکانس خاص است. رزولوشن فرکانسی توسط فرمول زیر محاسبه می‌شود:

Δf=fsN\Delta f = \frac{f_s}{N}

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

تصویر زیر به عنوان خروجی این کد ترسیم می‌شود.

تبدیل فوریه سریع
تبدیل فوریه سریع

همان طور که می‌بینیم، تبدیل FFT سیگنال به دست می‌آید و اطلاعاتی راجع به فرکانس‌های امواجی که در حوزه زمان بودند، در اختیار ما قرار می‌دهد. FFT یک مبادله بین اطلاعات در حوزه زمان و اطلاعات در حوزه فرکانس است. با گرفتن تبدیل FFT از یک سیگنال در حوزه زمان، تمام اطلاعات زمانی سیگنال را از دست داده و در عوض اطلاعات فرکانسی آن را به دست می‌آوریم. برای داشتن اطلاعات هم در حوزه زمان و هم در حوزه فرکانس در یک طیف، باید از طیف نگاره (Spectrogram) استفاده کنیم.

تصویر زیر FFT یک سیگنال موسیقی را نشان می‌دهد.

تبدیل فوریه سریع سیگنال موسیقی
تبدیل فوریه سریع سیگنال موسیقی

فرکانس این نمودار در مقیاس لگاریتمی ترسیم شده است. همان طور که در تصویر دیده می‌شود فرکانس غالب این سیگنال بین 101.510^{1.5} و 102.210^{2.2} هرتز (30 تا 158 هرتز) است. اگر مطابق مثال ابتدای متن، پنجره‌ای با فرکانس طبیعی 100 هرتز داشته باشیم، آن‌گاه می‌توان به این نکته پی برد که فرکانس طبیعی پنجره دقیقا در میان این بازه فرکانسی غالب سیگنال صوتی است و با بالا بردن صدای موسیقی، این امکان وجود دارد که پدیده رزونانس به وجود بیاید.

برنامه تبدیل فوریه سریع در جاوا

با استفاده از برنامه زیر می‌توان تبدیل فوریه سریع را در زبان برنامه‌نویسی جاوا محاسبه کرد.

تبدیل فوریه سریع در زبان متلب

با استفاده از برنامه زیر می‌توان تبدیل فوریه سریع را در زبان متلب محاسبه کرد.

تبدیل فوریه سریع در زبان C

با استفاده از برنامه زیر می‌توان تبدیل فوریه سریع را در زبان C محاسبه کرد.

تبدیل فوریه سریع در زبان #C

با استفاده از کد زیر، می‌توان تبدیل فوریه سریع را با زبان #C پیاده‌سازی کرد.

تبدیل فوریه سریع در زبان جاوا اسکریپت

جهت دسترسی به کدهای پیاده‌سازی تبدیل فوریه سریع یا fft در جاوا اسکریپت، روی این لینک کلیک کنید و فایل مربوطه حاوی کدهای جاوا اسکریپت را دریافت کنید.

اگر مطلب بالا برای شما مفید بوده و علاقه‌مند به موضوعات مرتبط با آن‌ هستید، پیشنهاد می‌کنیم آموزش‌های زیر را نیز ببینید:

^^

بر اساس رای ۵۳ نفر
آیا این مطلب برای شما مفید بود؟
اگر بازخوردی درباره این مطلب دارید یا پرسشی دارید که بدون پاسخ مانده است، آن را از طریق بخش نظرات مطرح کنید.
منابع:
ALL ABOUT CIRCUITSGitHubRitchie Vinknayuki
دانلود PDF مقاله
۸ دیدگاه برای «تبدیل فوریه سریع (fft) – به زبان ساده»

عالیه؛ واقعا خلاصه، مفید و تا حد زیادی کامل بود.
اگه امکانش هست بحث بیان FFT با ماتریس تنک که در بالا اشاره کردین ولی بهش نپرداختین رو به همراه بحث محاسبه انتگرال های فوریه با استفاده ازFFT (Computing Fourier Integrals Using the FFT) رو هم اضافه کنین. ممنون

ممنون. خدا خیرتان دهد

با درود
بسیار عالی و کارآمد بود.
سپاس فراوان

سلام…خداقوت بابت مطالب خوبتون.
من بامتلب اشنایی زیادی ندارم و برای پایان نامه باید از تبدیل فوریه استفاده کنم ، منتهی داده هایی از خارج محیط متلب رو وارد متلب کردم و نمیدونم چطور دیتاهای مربوطه رو وارد کدهای تبدیل فوریه که در بالا اورده شده ، کنم…ممنون میشم راهنمایی بفرمایید

عالی بود. ممنون.

ممنون از مطلب جامع و پر بهره تون…دست مریزاد گل کاشتید

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

با سلام و احترام
ممنون از مطالب خوبتون.
در طی مرور مطالب فوق، یک سوال برام پیش آمد در صورت امکان ممنون میشم راهنمایی بفرمایید.
در کدها بفرض مثال در متلب، ورودی های تابع راتوضیح می دید بایستی چی باشند.
function [spectrum, freq] = autofft(xs, ts, fftset
آیا xs ورودی های تابع در حوزه زمان
ts زمان نمونه برداری
fftset هم تابع مدنظره؟
با تشکر
حمید

نظر شما چیست؟

نشانی ایمیل شما منتشر نخواهد شد. بخش‌های موردنیاز علامت‌گذاری شده‌اند *