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

۸۷۷۹ بازدید
آخرین به‌روزرسانی: ۲۵ اردیبهشت ۱۴۰۲
زمان مطالعه: ۳۱ دقیقه
تبدیل فوریه سریع (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 تبدیل فوریه گسسته ترسیم شده است.

نمودار محاسبات مورد نیاز برای DFT با افزایش N
نمودار محاسبات مورد نیاز برای DFT با افزایش 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} $$ نقطه اولی در نمودار شکل بالا را می‌توان به صورت زیر نشان داد.

شکستن تبدیل فوریه گسسته $$ \tfrac{N}{2} $$ نقطه در نمودار قبل
شکستن تبدیل فوریه گسسته $$ \tfrac{N}{2} $$ نقطه در نمودار قبل

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

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

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

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

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

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

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

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

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

^^

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

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

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

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

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

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

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

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

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

نظر شما چیست؟

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