تبدیل فوریه گسسته — از صفر تا صد

در آموزشهای پیشین مجله فرادرس، با تبدیل فوریه و ویژگیهای آن آشنا شدیم. تبدیل فوریه گسسته یکی از ابزارهای قدرتمند در پردازش سیگنال دیجیتال است که با استفاده از آن میتوان طیف سیگنالی با طول متناهی را به دست آورد. در این آموزش، با تبدیل فوریه گسسته آشنا میشویم.
مواقع و شرایط زیادی وجود دارد که نیاز داریم محتوای فرکانسی یک سیگنال حوزه زمان را مشخص کنیم. برای مثال، ممکن است بخواهیم طیف خروجی یک نوسانساز LC را برای مشاهده مقدار نویز موجود در سیگنال سینوسی تولیدی تحلیل کنیم. این کار با استفاده از تبدیل فوریه گسسته (Discrete Fourier Transform) یا DFT میسر است. تبدیل فوریه گسسته، در کنار فیلترسازی دیجیتال، یکی از دو ابزار قوی پردازش سیگنال دیجیتال است.
تبدیل فوریه گسسته چیست؟
سیگنال پیوسته $$ x (t)$$ را در نظر بگیرید که در شکل ۱ نشان داده شده است. میخواهیم این سیگنال را تحلیل کنیم.

بدیهی است که یک کامپیوتر دیجیتال نمیتواند با یک سیگنال پیوسته کار کند و لازم است که تعدادی نمونه از $$ x (t)$$ بگیریم و به جای سیگنال اصلی، این نمونهها را تحلیل کنیم. علاوه بر این، علیرغم اینکه شکل ۱ تنها ۵ میلیثانیه از سیگنال را نشان میدهد، ممکن است $$ x(t)$$ ساعتها، سالها یا بیشتر از آن ادامه داشته باشد. از آنجایی که کامپیوتر دیجیتال توانایی فقط تعداد محدودی از نمونهها را دارد، باید از تقریب استفاده کرده و از تعداد محدودی نمونه بهره ببریم. بنابراین، عموماً، دنبالهای به طول محدود برای نمایش این سیگنال پیوسته آنالوگ به کار میرود که میتواند روی محور زمان به مثبت بینهایت میل کند. از سیگنال $$ x(t)$$ شکل ۱ با نرخ نمونهبرداری ۸۰۰۰ نمونه بر ثانیه، تعداد $$L=8$$ نمونه میگیریم. نتیجه این نمونهبرداری در شکل ۲ نشان داده شده است.

اگر محور زمان را به دوره نمونهبرداری $$T_s=\frac{1}{f_s}$$ بهنجار (نرمالیزه) کنیم، دنباله گسسته $$x(n)$$ را به دست خواهیم آورد که در جدول زیر ارائه شده است:
$$7$$ | $$6$$ | $$5$$ | $$4$$ | $$3$$ | $$2$$ | $$1$$ | $$0$$ | $$n$$ |
$$-0.8321$$ | $$-1.2165$$ | $$-0.5821$$ | $$0.2165$$ | $$0.5821$$ | $$0.7835$$ | $$0.8321$$ | $$0.2165$$ | $$x(n)$$ |
تحلیل فوریه چندین ابزار ریاضی را برای تعیین محتوای فرکانسی یک سیگنال حوزه زمان ارائه میکند، اما کدامیک از این ابزارها برای تحلیل $$ x (n)$$ مناسبتر است؟ تنها دو روش از خانواده تحلیل فوریه وجود دارد که مختص سیگنالهای گسستهاند: تبدیل فوریه زمانگسسته (Discrete-time Fourier Transform) یا DTFT و تبدیل فوریه گسسته (Discrete Fourier Transform) یا DFT.
تبدیل فوریه زمانگسسته (DTFT) دنباله ورودی $$ x(n)$$ به صورت زیر بیان میشود:
$$ \large X ( e ^ { j \omega } ) = \sum _ { n = – \infty } ^ { + \infty } x ( n ) e ^ { – j n \omega } \; \; \; \; \; ( 1 ) $$
معکوس DTFT نیز به شکل زیر است:
$$ \large x ( n ) = \frac { 1 } { 2 \pi } \int _ { – \pi } ^ { \pi } X ( e ^ { j \omega } ) e ^ { j n \omega } d \omega \; \; \; \; \; ( 2 ) $$
میتوانیم از معادله (۱) برای یافتن طیف سیگنال $$x(n)$$ استفاده کنیم که دامنه زمانی آن متناهی است. در معادله بالا، $$X(e^{j \omega})$$ تابعی پیوسته از $$ \omega $$ است. بنابراین، یک کمپیوتر دیجیتال نمیتواند مستقیماً از معادله (۱) برای تحلیل $$ x (n)$$ استفاده کند. البته میتوانیم از نمونههای $$X(e^{j \omega}) $$ برای یافتن تقریبی از طیف $$ x(n)$$ استفاده کنیم. ایده نمونهبرداری $$X(e^{j \omega}) $$ در نقاط فرکانسی با فاصله برابر، در واقع، اساس روش دوم فوریه، یعنی DFT است که در بالا به آن اشاره کردیم. توجه کنید که این نمونهبرداری در حوزه فرکانس انجام میشود ($$X(e^{j \omega})$$ تابعی از فرکانس است).
ریاضیات تبدیل فوریه گسسته (DFT)
این بخش را با بررسی مثالی که در بخش قبل بیان کردیم معرفی میکنیم. دنبالهای به طول $$L$$ از سیگنال $$ x(n)$$ داریم که مربوط سیگنال پیوسته آنالوگ $$x(t) $$ است. هدفی که به دنبال آن هستیم، یافتن مجموعهای از سینوسیها است که میتوان آنها را با هم جمع و $$ x(n)$$ را تولید کرد.
همانطور که در بخش قبلی گفتیم، DFT مبتنی بر نمونهبرداری در نقاط فرکانسی با فواصل برابر از DTFT است که در رابطه (۱) داده شد. $$ X(e^{j \omega}) $$ یک تابع دورهای یا متناوب از $$ \omega $$ با دوره تناوب $$2 \pi$$ است. اگر در هر دوره $$X(e^{j \omega})$$، تعداد $$N$$ نمونه بگیریم، فاصله بین نقاط نمونهبردای $$ \frac{ 2\pi }{N}$$ خواهد بود. بنابراین، فرکانس مجموعه سینوسیهایی که به دنبال یافتن آنها هستیم، به فرم $$\frac{2 \pi}{N} \times k $$ خواهد بود که در آن، $$k=0, 1, …, N-1$$. با استفاده از نماییهای مختلط مشابه با معادلات (۱) و (۲)، توابع پایه به صورت $$e^{j \frac{2 \pi}{N}kn}$$ خواهند بود. میخواهیم یک مجموع وزندار از این توابع را پیدا کنیم که سیگنال اصلی $$x(n)$$ را نتیجه دهند. در نتیجه، خواهیم داشت:
$$ \large x ( n ) = \sum _ { k = 0 } ^ { N – 1 } X’ ( k ) e ^ { j \frac { 2 \pi } { N } k n }, \; \; \; n=0, 1, \dots, L-1 \; \; \; \; \; ( 3 ) $$
که در آن، $$X'(k)$$ وزن نمایی مختلط $$ e^{j \frac{2 \pi}{N}kn}$$ است. معادله (۳)، معادل با مجموعه معادلات زیر است:
$$ \large
x ( 0 ) = X’ ( 0 ) + X’ ( 1 ) + X’ ( 2 ) + \dots + X’ ( N – 1 ) \\
\large x ( 1 ) = X’ ( 0 ) + X’ ( 1 ) e ^ { j \frac { 2 \pi }{ N } } + X’ ( 2 ) e ^ { j \frac { 2 \pi } { N } 2 } + \dots + X’ ( N – 1 ) e ^ { j \frac { 2 \pi } { N } ( N – 1 ) } \\
\large x ( 2 ) = X’ ( 0 ) + X’ ( 1 ) e ^ { j \frac { 2 \pi }{ N } 2 } + X’ ( 2 ) e ^ { j \frac { 2 \pi } { N } 2 \times 2 } + \dots + X’ ( N – 1 ) e ^ { j \frac { 2 \pi } { N } ( N – 1 ) \times 2 } \; \; \; \; \; \; \; \; (4) \\
\large \vdots \\
\large
x ( L – 1 ) = X’ ( 0 ) + X’ ( 1 ) e ^ { j \frac { 2 \pi } { N } ( L – 1 ) } + X’ ( 2 ) e^ { j \frac { 2 \pi } { N } 2 \times ( L – 1 ) } + \dots + X’ ( N – 1 ) e ^ { j \frac { 2 \pi } { N } ( N – 1 ) \times ( L – 1 ) } $$
برای یک $$L$$ و $$N$$ مشخص، مقادیر نماییهای مختلط معلوم بوده و با داشتن مقدار سیگنال حوزه زمان میتوانیم ضرایب $$X'(k) $$ را محاسبه کنیم. تعداد نمونههای حوزه زمان ($$L$$)، تعداد معادلات را در مجموعه معادلات و تعداد نمونههای حوزه فرکانس $$X(e^{j\omega})$$، تعداد مجهولات معادله (۴) را مشخص میکند. با قرار دادن $$L = N $$، تعداد $$N$$ معادله مستقل خطی برای یافتن $$N$$ مجهولِ $$X'(k) $$ خواهیم داشت.
با فرض $$L=N=8$$، میتوانیم فرم ماتریسی معادله (۴) را نوشته و با استفاده از متلب، مجهولات را به دست آوریم:
$$7$$ | $$6$$ | $$5$$ | $$4$$ | $$3$$ | $$2$$ | $$1$$ | $$0$$ | $$k$$ |
$$0.5j$$ | $$0.1083+0.0625j$$ | $$0$$ | $$0$$ | $$0$$ | $$0.1083-0.0625j$$ | $$-0.5j$$ | $$0$$ | $$X'(k)$$ |
اما تفسیر این تحلیل چیست؟ ما دنباله $$ x (n) $$ را به صورت مجموع نماییهای مختلط معادله (۳) نوشتیم. اکنون نتایج را در معادله مذکور جایگذاری میکنیم. با $$L=N=8$$ و حذف نماییهایی که ضریبشان صفر است، طبق جدول بالا میتوان نوشت:
$$ \large x ( n ) = X’ ( 1 ) e ^ { j \frac { 2 \pi }{ 8 } n } + X’ ( 2 ) e ^ { j \frac { 2 \pi } { 8 } 2 n } + X’ ( 6 ) e ^ { j \frac { 2 \pi } { 8 } 6 n } + X’ ( 7 ) e ^ { j \frac { 2 \pi } { 8 } 7 n } $$
از آنجایی که این توابع نمایی دورهای با $$N=8$$ هستند، میتوانیم این معادله را سادهتر کنیم. برای مثال، $$e^{j\frac{2\pi}{8}7n}$$ برابر با $$e^{j\frac{2\pi}{8}(8-1)n}=e^{j\frac{2\pi}{8}8n}e^{-j\frac{2\pi}{8}n}=e^{-j\frac{2\pi}{8}n} $$ است. بعد از سادهسازی و بازآرایی جملات معادله بالا، خواهیم داشت:
$$ \large x ( n ) = X’ ( 1 ) e ^ { j \frac { 2 \pi } { 8 }n } + X’ ( 7 ) e ^ {- j \frac { 2 \pi } { 8 } n } + X’ ( 2 ) e ^ { j \frac { 2 \pi }{ 8 } 2 n } + X’ ( 6 ) e ^ { – j \frac { 2 \pi } { 8 } 2 n } $$
در نهایت، رابطه زیر را به دست خواهیم آورد:
$$ \large \begin {align*}
x ( n ) & = \sin ( \frac { 2 \pi } { 8} n ) + 0 . 1 2 5 \sin ( \frac { 4 \pi } { 8 } n ) + 0 . 2 1 6 6 \cos ( \frac { 4 \pi }{ 8 } n ) \\ &= \sin ( \frac { 2 \pi } { 8 } n) + 0 . 2 5 \sin ( \frac { 4 \pi }{ 8 } n+ 3 \frac { \pi } { 3 } )
\end {align*} \; \; \; \; \; (5)$$
رابطه بالا نشان میدهد که میتوانیم $$x(n)$$ را با دو مؤلفه نشان دهیم: یکی در فرکانس نرمالیزه $$\frac{2\pi}{8}$$ با دامنه $$1$$ و دیگری در $$\frac{4\pi}{8} $$ با دامنه $$0.25$$. این فرکانسها نرمالیزه هستند؛ اگر فرکانس نمونهبرداری $$f_s=8000$$ نمونه بر ثانیه باشد، فرکانس این دو مؤلفه، به ترتیب، $$f_1=\frac{2\pi}{8}\frac{f_s}{2\pi}=1 kHz$$ و $$f_2=\frac{4\pi}{8}\frac{f_s}{2\pi}=2 kHz$$ خواهد بود.
اکنون خلاصهای را از آنچه که گفتیم بیان میکنیم. ما با یک سیگنال زمانپیوسته شروع کردیم و از تعداد محدودی نمونه برای تحلیل محتوای فرکانسی این سیگنال پیوسته استفاده کردیم. دیدیم که مسئله تجزیه این دنباله گسسته به حل مجموعهای از معادلات خطی میانجامد. در مثال سادهای که بررسی کردیم، معادله (۵) را نیز ساده کرده و یک معادله صریح برای $$ x (n)$$ به دست آوردیم.
این، نکتهای است که باید توجه ویژهای به آن داشته باشیم. تحلیل با دنبالهای به طول ۸ شروع میشود که در شکل ۳ نشان داده شده است. واضح است که مقدار $$x(n)$$ را برای $$ n=0, \dots, 7 $$ میدانیم، اما اطلاعی از مقدار آن در خارج از این محدوده نداریم.

تحلیل بالا، در واقع، جستوجو برای مجموعی از نماییهای مختلط است که میتوانند مقادیر $$ x(n)$$ را برای $$n=0, \dots, 7 $$ بازتولید کنند. معادله (۵) را در خارج از این بازه آزمایش میکنیم. برای جلوگیری از ابهام، به دنباله $$p(n)$$ بر میگردیم که در معادله (۵) داده شد. نمودار $$p(n)$$ که یک تابع دورهای با $$ N=8$$ است، در شکل ۴ نشان داده شده است. همانطور که در این شکل میبینیم، مقادیر $$ x (n)$$ اصلی هر ۸ نمونه یکبار تکرار میشوند. به عبارت دیگر، وقتی بتوانیم $$x(n)$$ را با استفاده از تعدادی نمایی مختلط نشان دهیم، نمایشِ به دست آمده متناوب بوده و $$x(n)$$ فقط در یک دوره تناوب برابر با $$p(n)$$ است.

توجه کنید که وقتی درباره تشکیل، تحلیل و اعمال DFT بحث میکنیم، منظورمان نمایشی از دنبالهای با طول متناهی با استفاده از یک دنباله متناوب است که در آن، مقادیرِ یک دوره تناوب از این دنباله دورهای برابر با مقادیر آن دنباله متناهی هستند.
همچنین دقت داشته باشید که رفتار دورهای $$p(n)$$ را میتوان با یادآوری این نکته درک کرد که در حال نمونهبرداری از $$X(e^{j\omega})$$ در حوزه فرکانس (تبدیل فوریه زمانگسسته $$x(n)$$) هستیم. در واقع، نمونهبرداری از $$ X(e^{j\omega})$$ در حوزه فرکانس، به نسخه برابر $$ x (n)$$ در حوزه زمان میانجامد. فرم تناوبی $$ x(n)$$ در شکل ۳ و $$p(n)$$ در شکل ۴ نشان داده شد.
استخراج معادلات DFT
روش بالا را که برای محاسبه طیف یک دنباله با طول متناهی بیان کردیم، ساده و قابل درک است و مبیّن رفتار تناوبی ذاتی نمایش DFT است. با وجود این، میتوانیم از مباحث بالا استفاده کرده و فرم بسته معادلات DFT را بدون نیاز به محاسبه معکوس یک ماتریس بزرگ به دست آوریم. برای این کار، تنها لازم است یک سیگنال تناوبی با $$N$$ نمونه از دنبالهای با طول متناهی بسازیم. پس از آن، با اعمال بسط سری فوریه زمانگسسته، میتوانیم نمایش حوزه فرکانس سیگنال تناوبی را به دست آوریم. ضرایب سری فوریه به دست آمده، جز برای یک عامل مقایسبندی، مشابه ضرایب DFT هستند.
فرض کنید دنباله متناهی که میخواهیم آن را تحلیل کنیم در شکل ۵ (الف) نشان داده شده است. برای محاسبه تبدیل فوریه گسسته Nنقطهای، لازم است از روی سیگنال $$x(n)$$ با دوره تناوب $$N$$، سیگنال $$P(n)$$ را بسازیم (شکل ۵ (ب)).

با توجه به رابطه $$p(n)=x(n)$$ به ازای $$n=0, 1, \dots, N-1$$ سری فوریه زمانگسسته این سیگنال متناوب به صورت زیر است:
$$ \large a _ k = \frac { 1 } { N } \sum _ { n = 0 } ^ { N – 1 } x ( n ) e ^ { – j \frac { 2 \pi } { N } k n } \; \; \; \; \; ( 6 ) $$
که در آن، $$N$$ دوره تناوب سیگنال است. سیگنال حوزه زمان را میتوان به صورت زیر به دست آورد:
$$ \large x ( n ) = \sum _ { k = 0 } ^ { N – 1 } a _ { k } e ^ { j \frac { 2 \pi } { N } k n } \; \; \; \; \; ( 7 ) $$
با ضرب ضرایب معادله (۶) در $$N$$، ضرایب تبدیل فوریه گسسته ($$X(n)$$) به دست میآید:
$$ \large X ( k ) = \sum _ { n = 0 } ^ { N – 1 } x ( n ) e ^ { -j \frac { 2 \pi } { N } k n } \; \; \; \; \; ( 8 ) $$
تبدیل فوریه گسسته معکوس به صورت زیر خواهد بود:
$$ \large x ( n ) = \frac { 1 } { N } \sum _ { k = 0 } ^ { N – 1 } X ( k ) e ^ { j \frac { 2 \pi } { N } k n } \; \; \; \; \; ( 9 ) $$
توجه کنید که وقتی سری فوریه زمانگسسته یک سیگنال متناوب باشد، ضرایب DFT یک دنباله متناهی در $$0 \leq k \leq N-1$$ هستند.
مثالها
در این بخش، دو مثال ساده را بررسی میکنیم.
مثال ۱
فرض کنید $$x_0=1$$، $$x_1=x_2= \cdots = x_{N-1} = 0 $$. تبدیل فوریه گسسته $$ x_n$$ به صورت زیر است:
$$ \large X _ k = \sum _ { n = 0 } ^ { N – 1 } x _ n e ^ { – 2 \pi i k n / N } = 1 . $$
بنابراین، $$x_n$$ به صورت زیر خواهد بود:
$$ \large x _ n = \frac 1 { N } \sum _ { k = 0 } ^ { N – 1 } e ^ { 2 \pi i k n / N } . $$
توجه کنید که وقتی $$n=0$$، $$ x _0 = \frac {1} {N} \cdot N = 1$$ و وقتی $$ n \neq 0$$، داریم:
$$ \large \begin {aligned} x _ n & = \frac 1 { N } \sum _ { k = 0 } ^ { N – 1 } \big ( e ^ { 2 \pi in / N } \big ) ^ k \\\\ & = \frac 1 { N } \frac { \big ( e ^ { 2 \pi in / N } \big ) ^ N – 1 } { e ^ { 2 \pi in / N } – 1 } \\\\ & = \frac 1 { N } \frac { e ^ { 2 \pi i n } – 1 } { e ^ { 2 \pi in / N } – 1 } \\\\ & = 0 . \end {aligned} $$
مثال ۲
تبدیل فوریه $$ f ( x _ 0 , x _ 1 , x _ 2 , x _ 3 ) = ( 0 , 1 , 0 , 0 ) $$ را به دست آورید.
حل: در این مثال، داریم:
$$ \large X _ k = \sum _ { n = 0 } ^ 3 x _ n e ^ { – 2 \pi i k n / 4 } = e ^ { – 2 \pi i k / 4 } . $$
بنابراین، میتوان گفت: $$X_0=1$$، $$X_1=-i$$، $$X_2 = -1$$ و $$X_3 = i$$. در نتیجه، پاسخ $$(1, -i , -1 , i )$$ است.
در نتیجه، $$x_n$$ به صورت زیر خواهد بود:
$$ \large x _ n = 1 – i e ^ { 2 \pi i n / 4 } – e ^ { 4 \pi i n / 4 } + i e ^ { 6 \pi i n / 4 } . $$
با تبدیل ضرایب مختلط به نماییهای مختلط، داریم:
$$ \large \begin {aligned} x _ n & = 1 + e ^ { 6 \pi i / 4 } e ^ { 2 \pi i n / 4 } + e ^ { 4 \pi i / 4 } e ^ { 4 \pi i n / 4 } + e ^ { 2 \pi i / 4 } e ^ { 6 \pi i n / 4 } \\ & = 1 + e ^ { 2 \pi i ( n + 3 ) / 4 } + e ^ { 2 \pi i ( 2 n + 2 ) / 4 } + e ^ { 2 \pi i ( 3 n + 1 ) / 4 } . \end {aligned} $$
این مثال، موردی از جابهجایی فاز است که اتفاق میافتد. با در نظر گرفتن بخشهای حقیقی هر دو سمت، مجموعی از موجهای کسینوسی خواهیم داشت:
$$ \large x _ n = 1 + \cos ( 2 \pi n / 4 + 3 \pi / 2 ) + \cos ( 4 \pi n / 4 + \pi ) + \cos ( 6 \pi n / 4 + \pi / 2 ) $$
که در آن، افزوده شدن $$3\pi/2$$، $$\pi$$ و $$\pi/2$$ موجب جابهجایی شکل موجها به ترتیب، به اندازه $$270^\circ$$، $$180^ \circ$$ و $$ 90 ^ \circ $$ خواهد شد.
تعامد و تبدیل معکوس
چرا فرمول تبدیل معکوس درست است؟ $$X_k$$ را در فرمول $$ x _ n$$ جایگذاری میکنیم:
$$ \large \begin {aligned} \sum _ { k = 0 } ^ { N – 1 } X _ k e ^ { – 2 \pi i k n / N } & = \frac 1 { N } \sum _ { k = 0 } ^ { N – 1 } \sum _ { m = 0 } ^ { N – 1 } x _ m e ^ { 2 \pi i k m / N } e ^ { – 2 \pi i k n / N } \\ & = \frac 1 { N } \sum _ { k= 0 } ^ { N – 1 } \sum _ { m = 0 } ^ { N – 1 } x _ m e ^ { 2 \pi i k ( m – n ) / N } \\ & = \frac 1 { N } \sum _ { m = 0 } ^ { N – 1 } x _ m \sum _ { k = 0 } ^ { N – 1 } e ^ { 2 \pi i k ( m – n ) / N } . \end {aligned} $$
اگر $$ m \neq n$$، مجموع داخلی طبق فرمول سری هندسی برابر با صفر است. وقتی $$ m = n$$، مجموع داخلی برابر با $$ N $$ خواهد بود. بنابراین، مجموع کلی برابر است با $$ \frac {1}{N} \cdot x _ n \cdot N = x _ n $$.
یک راه دیگر برای بررسی این استدلال، آن است که این امر نتیجه تعامد نسبت به ضرب نقطهای مختلط است. بردارهای مختلط زیر را در نظر بگیرید:
$$ \large v _ k = \left ( 1 , \omega _ N ^ k , \omega _ N ^{ 2 k } , \ldots , \omega _ N ^ { k ( N – 1 ) } \right ) , $$
که در آنها، $$ \omega _ N = e ^ {2 \pi i /N}$$ است. با استدلالی مشابه آنچه در بالا گفتیم، این بردارهای مختلط نسبت به ضرب نقطهای مختلط متعامد هستند: $$ v _k \cdot v _ l = N$$ اگر $$ k = l $$ و در غیر این صورت، برابر با صفر است.
فرمول DFT برای $$ X_k$$ به سادگی به صورت $$ X _ k = x \cdot v _ k $$ بیان میشود، که در آن، $$x$$ بردار $$ ( x _ 0, x _ 1, \cdots , x _ {N-1})$$ است. فرمول معکوس، $$ x $$ را با نوشتن $$x$$ با استفاده از فرمول استاندارد برای بیان هر بردار به عنوان ترکیب خطی بردارهای پایه متعامد بر حسب $$X$$ بر میگرداند:
$$ \large x = \sum _ { k = 0 } ^ { N – 1 } \frac { x \cdot v _ k }{ v _ k \cdot v _ k } v _ k = \frac 1 { N } \sum _ { k = 0 } ^ { N – 1 } X _ k v _ k . $$
پیادهسازی تبدیل فوریه گسسته
در این بخش، مراحل پیادهسازی الگوریتم تبدیل فوریه گسسته را بیان کرده و آن را در زبانهای برنامهنویسی مختلف ارائه میکنیم.
در حالت کلی میتوان گفت که DFT تابعی است که برداری با $$n$$ عدد مختلط را به برداری با $$n$$ عدد مختلط دیگر مینگارد. با در نظر گرفتن مبداء $$0$$، فرض میکنیم $$ x (t)$$ نشان دهنده $$t$$اُمین درایه بردار ورودی و $$ X(k)$$ نشان دهنده $$k$$اُمین درایه بردار خروجی باشد. تبدیل فوریه گسسته پایه با فرمول زیر بیان میشود:
$$ \large X ( k ) = \displaystyle \sum _ { t = 0 } ^ { n – 1 } x ( t ) e ^ { – 2 \pi i t k / n } $$
فرمول بالا مبیّن این موضوع است که بردار $$x$$ سطح یا اندازه سیگنال را در نقاط زمانی مختلف نشان میدهد و بردار $$X$$ مجموع اندازه را در فرکانسهای مختلف به نمایش میگذارد. آنچه این فرمول میگوید، این است که اندازه سیگنال در فرکانس $$k$$ برابر است با مجموع اندازههای سیگنال در هر زمان $$t$$ ضرب در یک نمایی مختلط.
با داشتن دانش کمی درباه نماد مجموع، اعداد مختلط و برنامهنویسی کامپیوتری میتوانیم به سادگی برنامه تبدیل فوریه گسسته را بنویسیم. در ادامه، این فرایند را گام به گام بررسی میکنیم. برای نمونه، این مراحل را در قالب زبان برنامهنویسی جاوا بیان خواهیم کرد.
ساختار برنامه
همانطور که گفتیم، سری فوریه گسسته یک بردار با $$n$$ عدد مختلط را به عنوان ورودی میگیرد و بردار خروجی را که شامل $$n$$ عدد مختلط است محاسبه میکند. از آنجایی که جاوا عدد مختلط ندارد، باید خودمان عدد مختلط را با استفاده از دو عدد حقیقی تعریف کنیم. بردار، دنبالهای از اعداد است که میتوان آن را با یک آرایه نمایش داد. به جای بازگرداندن آرایههای خروجی، آنها را به عنوان آرگومان در مرجع قرار خواهیم داد. ساختار آن به صورت زیر است:
void dft(double[] inreal , double[] inimag, double[] outreal, double[] outimag) { // Assume all 4 arrays have the same length int n = inreal.length; // Incomplete: Method body }
در ادامه حلقه بیرونی را برای اختصاص دادن یک مقدار به هر درایه خروجی مینویسیم:
void dft(double[] inreal , double[] inimag, double[] outreal, double[] outimag) { int n = inreal.length; for (int k = 0; k < n; k++) { // For each output element outreal[k] = ?; // Incomplete outimag[k] = ?; // Incomplete } }
مجموع
نوشتن مجموع یا همان سری شاید در ابتدا کمی پیچیده به نظر برسد، اما در واقع برعکس است. فرم کلی یک مجموع به صورت زیر ساده میشود:
$$ \large \displaystyle \sum _ { j = a } ^ { b } f ( j ) = f ( a ) + f ( a + 1 ) + \cdots + f ( b – 1 ) + f ( b ) $$
در حقیقت، در سری بالا مقادیر $$j$$ را یک به یک جایگذاری میکنیم. به سادگی کد سری بالا را میتوان به صورت زیر نوشت:
double sum = 0; for (int j = a; j <= b; j++) { sum += f(j); } // The value of 'sum' now has the desired answer
ریاضیات عدد مختلط
جمغ دو عدد مختلط به سادگی انجام میشود:
$$ \large ( a + b i ) + ( c + d i ) = ( a + c ) + ( b +d ) i $$
ضرب دو عدد مختلط کمی سختتر است؛ با استفاده از قانون توزیع و اتحاد $$ i ^ 2 = – 1 $$، داریم:
$$ \large ( a + b i ) ( c + d i ) = a c + a d \, i + b c \, i – b d = ( a c – b d ) + ( a d + b c ) i $$
طبق رابطه اویلر، برای هر عدد حقیقی $$ x$$، میتوان تساوی $$ e^{xi} = \cos x + i \sin x $$ را نوشت. علاوه بر این، کسینوس یک تابع زوج بوده و رابطه $$ \cos(-x) = \cos x $$ برقرار است. همچنین، از آنجایی که سینوس یک تابع فرد است، میتوان تساوی $$\sin(-x) = -(\sin x) $$ را نوشت. بنابراین، خواهیم داشت:
$$ \large \begin {align*} e ^ { – 2 \pi i t k / n } = e ^ { ( – 2 \pi t k / n ) i } & = \cos \left ( – 2 \pi \frac { t k } { n } \right ) + i \sin \left ( -2 \pi \frac { t k } { n } \right ) \\ &= \cos \left ( 2 \pi \frac { t k } { n } \right ) – i \sin \left ( 2 \pi \frac { t k } { n } \right ) \end {align*} $$
نماد $$\text{Re}(x)$$ را به عنوان بخش حقیقی $$x$$ و $$ \text{Im}(x) $$ را به عنوان بخش موهومی آن در نظر میگیریم. طبق تعریف، $$x = \text{Re}(x) + i \, \text{Im}(x) $$. بنابراین:
$$ \large x ( t ) \, e ^ { – 2 \pi i t k / n } = \left [ \text {Re} ( x ( t ) ) + i \, \text {Im} ( x ( t ) ) \right ] \left [ \cos \left ( 2 \pi \frac { t k } { n } \right ) – i \sin \left ( 2 \pi \frac { t k } { n } \right ) \right ] $$
با گسترش ضرب مختلط، داریم:
$$ \large \begin {align*} x ( t ) \, e ^ { – 2 \pi i t k / n } & = \left [ \text {Re} ( x ( t ) ) \cos \left ( 2 \pi \frac { t k } { n } \right ) + \text {Im}( x ( t ) ) \sin \left ( 2 \pi \frac { t k } { n } \right ) \right ] \\ & \;\;\;\; + i \left [ -\text {Re} ( x ( t ) ) \sin \left ( 2 \pi \frac { t k } { n } \right ) + \text {Im} ( x ( t ) ) \cos \left ( 2 \pi \frac { t k } { n } \right ) \right ] \end {align*} $$
بنابراین، هر عبارت در مجموع کد زیر را برای بخشهای حقیقی و موهومی خواهد داشت:
double angle = 2 * Math.PI * t * k / n; double real = inreal[t] * Math.cos(angle) + inimag[t] * Math.sin(angle); double imag = -inreal[t] * Math.sin(angle) + inimag[t] * Math.cos(angle);
ترکیب همه بخشها با یکدیگر
اگر همه بخشها را که گفتیم، با هم ترکیب کنیم، برنامه زیر را برای محاسبه سری فوریه گسسته خواهیم داشت:
static void dft(double[] inreal , double[] inimag, double[] outreal, double[] outimag) { int n = inreal.length; for (int k = 0; k < n; k++) { // For each output element double sumreal = 0; double sumimag = 0; for (int t = 0; t < n; t++) { // For each input element double angle = 2 * Math.PI * t * k / n; sumreal += inreal[t] * Math.cos(angle) + inimag[t] * Math.sin(angle); sumimag += -inreal[t] * Math.sin(angle) + inimag[t] * Math.cos(angle); } outreal[k] = sumreal; outimag[k] = sumimag; } }
در ادامه، کد DFT را برای چند زبان برنامهنویسی مختلف ارائه میکنیم. قبل از بررسی این کدها، نکات زیر را در نظر داشته باشید:
- زبانهای برنامهنویسی پایتون، سی، سیپلاسپلاس، سیشارپ، و متلب، اعداد مختلط را به صورت توکار دارند. این ویژگی کارمان را برای نوشتن برنامه DFT بسیار آسانتر میکند.
- هریک از برنامهها، نامگذاری قراردادی، سبک قالببندی و اصطلاحات کدنویسی مربوط به خود را دارند و لزوماً شبیه آن چیزی نیستند که برای جاوا بیان کردیم.
- تبدیل فوریه گسستهای که در این بخش بیان کردیم، برای سادگی، شامل نرمالیزه کردن/مقیاسبندی نیست. گزینههای متعددی برای چگونگی مقیاسبندی DFT وجود دارد. مثلاً مقیاسبندی فقط تبدیل مستقیم با $$1/n$$، مقیاسبندی تبدیلهای مستقیم و معکوس با $$1/\sqrt{n}$$، مقیاسبندی جداول مثلثاتی از پیش محاسبه شده و… .
برنامه تبدیل فوریه گسسته در جاوا
/* * Discrete Fourier transform (Java) * by Project Nayuki, 2017. Public domain. * https://www.nayuki.io/page/how-to-implement-the-discrete-fourier-transform */ public final class Dft { /* * Computes the discrete Fourier transform (DFT) of the given complex vector. * All the array arguments must be non-null and have the same length. */ public static void computeDft(double[] inreal, double[] inimag, double[] outreal, double[] outimag) { int n = inreal.length; for (int k = 0; k < n; k++) { // For each output element double sumreal = 0; double sumimag = 0; for (int t = 0; t < n; t++) { // For each input element double angle = 2 * Math.PI * t * k / n; sumreal += inreal[t] * Math.cos(angle) + inimag[t] * Math.sin(angle); sumimag += -inreal[t] * Math.sin(angle) + inimag[t] * Math.cos(angle); } outreal[k] = sumreal; outimag[k] = sumimag; } } }
برنامه تبدیل فوریه گسسته در C
/* * Discrete Fourier transform (C) * by Project Nayuki, 2017. Public domain. * https://www.nayuki.io/page/how-to-implement-the-discrete-fourier-transform */ /* * Computes the discrete Fourier transform (DFT) of the given complex vector. * All the array arguments must be non-NULL and have a length equal to n. */ #include <complex.h> #include <math.h> void compute_dft_complex(const double complex input[], double complex output[], int n) { for (int k = 0; k < n; k++) { // For each output element complex double sum = 0.0; for (int t = 0; t < n; t++) { // For each input element double angle = 2 * M_PI * t * k / n; sum += input[t] * cexp(-angle * I); } output[k] = sum; } } /* * (Alternate implementation using only real numbers.) * Computes the discrete Fourier transform (DFT) of the given complex vector. * All the array arguments must be non-NULL and have a length equal to n. */ #include <math.h> void compute_dft_real_pair(const double inreal[], const double inimag[], double outreal[], double outimag[], int n) { for (int k = 0; k < n; k++) { // For each output element double sumreal = 0; double sumimag = 0; for (int t = 0; t < n; t++) { // For each input element double angle = 2 * M_PI * t * k / n; sumreal += inreal[t] * cos(angle) + inimag[t] * sin(angle); sumimag += -inreal[t] * sin(angle) + inimag[t] * cos(angle); } outreal[k] = sumreal; outimag[k] = sumimag; } }
برنامه تبدیل فوریه گسسته در C++
/* * Discrete Fourier transform (C++) * by Project Nayuki, 2017. Public domain. * https://www.nayuki.io/page/how-to-implement-the-discrete-fourier-transform */ // Shared definitions #include <cmath> #include <vector> using std::size_t; using std::vector; /* * Computes the discrete Fourier transform (DFT) of the given complex vector. * All the array arguments must have the same length. */ #include <complex> using std::complex; using std::exp; vector<complex<double> > computeDft(const vector<complex<double> > &input) { vector<complex<double> > output; size_t n = input.size(); for (size_t k = 0; k < n; k++) { // For each output element complex<double> sum(0.0, 0.0); for (size_t t = 0; t < n; t++) { // For each input element double angle = 2 * M_PI * t * k / n; sum += input[t] * exp(-angle); } output.push_back(sum); } return output; } /* * (Alternate implementation using only real numbers.) * Computes the discrete Fourier transform (DFT) of the given complex vector. * All the array arguments must have the same length. */ using std::cos; using std::sin; void computeDft(const vector<double> &inreal, const vector<double> &inimag, vector<double> &outreal, vector<double> &outimag) { size_t n = inreal.size(); for (size_t k = 0; k < n; k++) { // For each output element double sumreal = 0; double sumimag = 0; for (size_t t = 0; t < n; t++) { // For each input element double angle = 2 * M_PI * t * k / n; sumreal += inreal[t] * cos(angle) + inimag[t] * sin(angle); sumimag += -inreal[t] * sin(angle) + inimag[t] * cos(angle); } outreal[k] = sumreal; outimag[k] = sumimag; } }
برنامه تبدیل فوریه گسسته در C#
/* * Discrete Fourier transform (C#) * by Project Nayuki, 2017. Public domain. * https://www.nayuki.io/page/how-to-implement-the-discrete-fourier-transform */ using System; using System.Numerics; // For Complex public sealed class Dft { /* * Computes the discrete Fourier transform (DFT) of the given complex vector. */ public static Complex[] computeDft(Complex[] input) { int n = input.Length; Complex[] output = new Complex[n]; for (int k = 0; k < n; k++) { // For each output element Complex sum = 0; for (int t = 0; t < n; t++) { // For each input element double angle = 2 * Math.PI * t * k / n; sum += input[t] * Complex.Exp(new Complex(0, -angle)); } output[k] = sum; } return output; } /* * Computes the discrete Fourier transform (DFT) of the given complex vector. * All the array arguments must be non-null and have the same length. */ public static void computeDft(double[] inreal, double[] inimag, double[] outreal, double[] outimag) { int n = inreal.Length; for (int k = 0; k < n; k++) { // For each output element double sumreal = 0; double sumimag = 0; for (int t = 0; t < n; t++) { // For each input element double angle = 2 * Math.PI * t * k / n; sumreal += inreal[t] * Math.Cos(angle) + inimag[t] * Math.Sin(angle); sumimag += -inreal[t] * Math.Sin(angle) + inimag[t] * Math.Cos(angle); } outreal[k] = sumreal; outimag[k] = sumimag; } } }
برنامه تبدیل فوریه گسسته در JavaScript
/* * Discrete Fourier transform (JavaScript) * by Project Nayuki, 2018. Public domain. * https://www.nayuki.io/page/how-to-implement-the-discrete-fourier-transform */ "use strict"; /* * Computes the discrete Fourier transform (DFT) of the given complex vector. * 'inreal' and 'inimag' are each an array of n floating-point numbers. * Returns an array of two arrays - outreal and outimag, each of length n. */ function computeDft(inreal, inimag) { var n = inreal.length; var outreal = new Array(n); var outimag = new Array(n); for (var k = 0; k < n; k++) { // For each output element var sumreal = 0; var sumimag = 0; for (var t = 0; t < n; t++) { // For each input element var angle = 2 * Math.PI * t * k / n; sumreal += inreal[t] * Math.cos(angle) + inimag[t] * Math.sin(angle); sumimag += -inreal[t] * Math.sin(angle) + inimag[t] * Math.cos(angle); } outreal[k] = sumreal; outimag[k] = sumimag; } return [outreal, outimag]; }
برنامه تبدیل فوریه گسسته در TypeScript
/* * Discrete Fourier transform (TypeScript) * by Project Nayuki, 2018. Public domain. * https://www.nayuki.io/page/how-to-implement-the-discrete-fourier-transform */ "use strict"; /* * Computes the discrete Fourier transform (DFT) of the given complex vector. * 'inreal' and 'inimag' are each an array of n floating-point numbers. * Returns an array of two arrays - outreal and outimag, each of length n. */ function computeDft(inreal: Array<number>, inimag: Array<number>): [Array<number>,Array<number>] { const n: number = inreal.length; let outreal: Array<number> = new Array(n); let outimag: Array<number> = new Array(n); for (let k = 0; k < n; k++) { // For each output element let sumreal: number = 0; let sumimag: number = 0; for (let t = 0; t < n; t++) { // For each input element const angle: number = 2 * Math.PI * t * k / n; sumreal += inreal[t] * Math.cos(angle) + inimag[t] * Math.sin(angle); sumimag += -inreal[t] * Math.sin(angle) + inimag[t] * Math.cos(angle); } outreal[k] = sumreal; outimag[k] = sumimag; } return [outreal, outimag]; }
برنامه تبدیل فوریه گسسته در Python
# # Discrete Fourier transform (Python) # by Project Nayuki, 2017. Public domain. # https://www.nayuki.io/page/how-to-implement-the-discrete-fourier-transform # # # Computes the discrete Fourier transform (DFT) of the given complex vector. # 'input' is a sequence of numbers (integer, float, or complex). # Returns a list of complex numbers as output, having the same length. # import cmath def compute_dft_complex(input): n = len(input) output = [] for k in range(n): # For each output element s = complex(0) for t in range(n): # For each input element angle = 2j * cmath.pi * t * k / n s += input[t] * cmath.exp(-angle) output.append(s) return output # # (Alternate implementation using only real numbers.) # Computes the discrete Fourier transform (DFT) of the given complex vector. # 'inreal' and 'inimag' are each a sequence of n floating-point numbers. # Returns a tuple of two lists of floats - outreal and outimag, each of length n. # import math def compute_dft_real_pair(inreal, inimag): assert len(inreal) == len(inimag) n = len(inreal) outreal = [] outimag = [] for k in range(n): # For each output element sumreal = 0.0 sumimag = 0.0 for t in range(n): # For each input element angle = 2 * math.pi * t * k / n sumreal += inreal[t] * math.cos(angle) + inimag[t] * math.sin(angle) sumimag += -inreal[t] * math.sin(angle) + inimag[t] * math.cos(angle) outreal.append(sumreal) outimag.append(sumimag) return (outreal, outimag)
برنامه تبدیل فوریه گسسته در Rust
/* * Discrete Fourier transform (Rust) * by Project Nayuki, 2017. Public domain. * https://www.nayuki.io/page/how-to-implement-the-discrete-fourier-transform */ fn compute_dft(inreal: &[f64], inimag: &[f64], outreal: &mut [f64], outimag: &mut [f64]) { let n: usize = inreal.len(); for k in 0 .. n { // For each output element let mut sumreal: f64 = 0.0; let mut sumimag: f64 = 0.0; for t in 0 .. n { // For each input element let angle: f64 = 2.0 * std::f64::consts::PI * (t as f64) * (k as f64) / (n as f64); sumreal += inreal[t] * angle.cos() + inimag[t] * angle.sin(); sumimag += -inreal[t] * angle.sin() + inimag[t] * angle.cos(); } outreal[k] = sumreal; outimag[k] = sumimag; } }
برنامه تبدیل فوریه گسسته در MATLAB
% % Discrete Fourier transform (MATLAB) % by Project Nayuki, 2017. Public domain. % https://www.nayuki.io/page/how-to-implement-the-discrete-fourier-transform % % % Computes the discrete Fourier transform (DFT) of the given complex vector. % 'input' can be a row vector or a column vector. The returned output has the same dimensions. % function output = compute_dft_scalarized(input) assert(isvector(input)); n = numel(input); output = NaN(size(input)); for k = 0 : n - 1 % For each output element s = 0; for t = 0 : n - 1 % For each input element s = s + input(t + 1) * exp(-2i * pi * t * k / n); end output(k + 1) = s; end end % % (Alternate implementation using matrix arithmetic.) % Computes the discrete Fourier transform (DFT) of the given complex vector. % 'input' can be a row vector or a column vector. The returned output has the same dimensions. % function output = compute_dft_vectorized(input) assert(isvector(input)); n = numel(input); matrix = exp(-2i * pi / n * (0 : n-1)' * (0 : n-1)); if size(input, 1) == 1 % Row vector output = input * matrix; elseif size(input, 2) == 1 % Column vector output = matrix * input; end end
جمعبندی
- تبدیل فوریه گسسته یکی از قویترین ابزارهای پردازش سیگنال دیجیتال است که به ما این امکان را میدهد طیف سیگنال متناهی $$ x(n)$$ را پیدا کنیم.
- اساساً، محاسبه DFT معادل با حل مجموعهای از معادلات خطی است.
- تبدیل فوریه گسسته با استفاده از یک دنباله متناوب نمایشی از یک دنباله متناهی در اختیار ما قرار میدهد. یک دوره تناوب دنباله متناوب مشابه با دنباله متناهی است. در نتیجه، میتوانیم از سری فوریه زمانگسسته برای استخراج معادلات DFT استفاده کنیم.
اگر مطلب بالا برای شما مفید بوده و علاقهمند به موضوعات مرتبط با آن هستید، پیشنهاد میکنیم آموزشهای زیر را نیز ببینید:
- مجموعه آموزشهای پردازش سیگنال
- آموزش پردازش سیگنالهای دیجیتال با متلب
- مجموعه آموزشهای مهندسی الکترونیک
- آموزش تجزیه و تحلیل سیگنال ها و سیستم ها
- تبدیل فوریه سریع (fft) — به زبان ساده
- تابع تبدیل پالسی — از صفر تا صد
- سری فوریه مختلط — به زبان ساده
^^
با سلام
بسیار عالی
فقط یک باید بگم در ابتدای بحث، تعداد نمونه برداری شما 1600 نقطه بر ثانیه است یعنی 1600 هرتز است. که اشتباها تایپ شده 8000 هرتز.
زیرا کل شکل 1، 5 میلی ثانیه است و 8 نمونه گرفته شده
8 تقسیم بر 5 میلی ثانیه میشه.. 1600 هرتز
با سلام؛
متن بازبینی شد و صحیح است.
با تشکر از همراهی شما با مجله فرادرس
چقدر روان مثال زدید و توضیح دادید… یکی از مباحث مهم در پردازش سیگنال نحوه پیاده سازی تبدیل فوریه هست که به ساده ترین شکل اون رو انجام دادید. سپاسگزارم
سلام.
خوشحالیم که این آموزش برایتان مفید بوده است.
شاد و پیروز باشید.
سلام
روز بخیر
عااالی بود کاملا متوجه شدم.
ممنون