آموزش تبدیل فوریه، DFT و FFT در یادگیری ماشین

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

این ابزار ریاضی پایه‌ی بسیاری از پیشرفت‌های یادگیری ماشین امروزی است؛ از اسپکتروگرام‌ها در مدل‌های تبدیل گفتار به متن مانند Whisper گرفته تا لایه‌های کانولوشن سریع در شبکه‌های بینایی و حتی کدگذاری مکانی سینوسی در ترنسفورمرها. در این درس، تبدیل فوریه‌ی گسسته (DFT) و الگوریتم سریع آن (FFT) را از صفر پیاده‌سازی می‌کنید تا درک شهودی و عمیقی از این مفاهیم به دست آورید.

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

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

اهداف یادگیری

• پیاده‌سازی DFT از صفر و اعتبارسنجی آن در مقایسه با الگوریتم Cooley-Tukey FFT با مرتبه‌ی $O(N \log N)$
• تفسیر ضرایب فرکانسی: استخراج دامنه، فاز و طیف توان از یک سیگنال
• به‌کارگیری قضیه‌ی کانولوشن برای انجام کانولوشن از طریق ضرب در دامنه‌ی فرکانس با FFT
• پیوند میان تجزیه‌ی فرکانسی فوریه با کدگذاری‌های مکانی ترنسفورمر و لایه‌های کانولوشن در CNN

صورت مسئله

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

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

تبدیل فوریه داده‌ها را از دامنه‌ی زمان به دامنه‌ی فرکانس منتقل می‌کند. این تبدیل، سیگنال را به امواج سینوسی با فرکانس‌های گوناگون تجزیه می‌کند. هر موج سینوسی یک دامنه (میزان شدت) و یک فاز (نقطه‌ی شروع) دارد؛ تبدیل فوریه هر دوی این اطلاعات را به شما می‌دهد.

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

مفاهیم پایه

تعریف DFT

به ازای $N$ نمونه به صورت $x[0], x[1], \ldots, x[N-1]$، تبدیل فوریه‌ی گسسته (DFT) تعداد $N$ ضریب فرکانسی $X[0], X[1], \ldots, X[N-1]$ تولید می‌کند:

$$X[k] = \sum_{n=0}^{N-1} x[n] \, e^{-2\pi i k n/N} \qquad \text{for } k = 0, 1, \ldots, N-1$$

هر $X[k]$ یک عدد مختلط است. اندازه‌ی آن یعنی $\lvert X[k] \rvert$، متناسب با دامنه‌ی فرکانس $k$ است. بدون نرمال‌سازی، یک سینوسی حقیقی با دامنه‌ی $A$ که دقیقاً روی خانه‌ی فرکانسی (Bin) شماره‌ی $k$ بیفتد مقدار $\lvert X[k] \rvert = AN/2$ را به دست می‌دهد؛ بنابراین برای خواندن دامنه باید آن را بر $N/2$ تقسیم کنید. فاز آن یعنی $\angle X[k]$، آفست فاز همان فرکانس را به شما می‌گوید.

نکته‌ی کلیدی: عبارت $e^{-2\pi i k n/N}$ یک فازور در حال چرخش با فرکانس $k$ است. تبدیل DFT در واقع همبستگی میان سیگنال و تک‌تکِ $N$ فرکانس با فواصل یکسان را محاسبه می‌کند. اگر سیگنال در فرکانس $k$ انرژی داشته باشد، این همبستگی بزرگ خواهد بود؛ در غیر این صورت مقداری نزدیک به صفر خواهد داشت.

معنای هر ضریب

$X[0]$: مؤلفه‌ی DC. این ضریب برابر با مجموع همه‌ی نمونه‌هاست و با میانگین تناسب دارد. این مؤلفه نشان‌دهنده‌ی آفست ثابت (فرکانس صفر) سیگنال است.

$$X[0] = \sum_{n=0}^{N-1} x[n] \, e^0 = \text{sum of all samples}$$

$X[k]$ به ازای $1 \le k \le N/2$: فرکانس‌های مثبت. ضریب $X[k]$ نمایانگر فرکانسِ $k$ چرخه در هر $N$ نمونه است. مقادیر بالاتر $k$ یعنی فرکانس بالاتر (نوسان سریع‌تر).

$X[N/2]$: فرکانس نایکویست (Nyquist). بالاترین فرکانسی است که می‌توانید با $N$ نمونه نمایش دهید. بالاتر از این مقدار دچار الیاسینگ (Aliasing) می‌شوید؛ پدیده‌ای که در آن فرکانس‌های بالا خود را به جای فرکانس‌های پایین جا می‌زنند.

$X[k]$ به ازای $N/2 < k < N$: فرکانس‌های منفی. برای سیگنال‌های حقیقی‌مقدار، رابطه‌ی $X[N-k] = \overline{X[k]}$ برقرار است. فرکانس‌های منفی بازتاب آینه‌ای فرکانس‌های مثبت هستند؛ به همین دلیل تمام اطلاعات مفید در همان $N/2 + 1$ ضریب اول قرار دارد.

تبدیل معکوس DFT

تبدیل معکوس DFT سیگنال اصلی را از روی ضرایب فرکانسی‌اش بازسازی می‌کند:

$$x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] \, e^{2\pi i k n/N} \qquad \text{for } n = 0, 1, \ldots, N-1$$

تنها تفاوت‌های آن با DFT مستقیم: علامت توان مثبت است (نه منفی) و یک فاکتور نرمال‌سازی $1/N$ نیز وجود دارد.

تبدیل معکوس DFT سیگنال را کاملاً و بدون نقص بازسازی می‌کند و هیچ داده‌ای از دست نمی‌رود. می‌توانید بدون هیچ خطایی از دامنه‌ی زمان به دامنه‌ی فرکانس بروید و برگردید. تبدیل DFT در حقیقت یک تغییر پایه (Change of basis) است؛ یعنی همان اطلاعات را در یک دستگاه مختصات دیگر بازنمایی می‌کند.

الگوریتم FFT: سرعت بخشیدن به محاسبات

محاسبه‌ی DFT بر اساس تعریف بالا از مرتبه‌ی $O(N^2)$ است: به ازای هر یک از $N$ ضریب خروجی، باید روی $N$ نمونه‌ی ورودی جمع بزنید. برای مقدار $N$ = 1 میلیون، این یعنی $10^{12}$ عملیات.

الگوریتم تبدیل فوریه‌ی سریع (FFT) همین خروجی را در زمان $O(N \log N)$ به دست می‌آورد. برای $N$ = 1 میلیون، این یعنی حدود 20 میلیون عملیات به جای یک تریلیون. همین جهش است که تحلیل فرکانس را به ابزاری کاربردی تبدیل می‌کند.

الگوریتم کولی-توکی (Cooley-Tukey)، به عنوان رایج‌ترین روش FFT، با رویکرد تقسیم و حل (Divide and conquer) کار می‌کند:

1. سیگنال را به دو بخشِ نمونه‌های با اندیس زوج و فرد تقسیم می‌کند.
2. تبدیل DFT هر نیمه را به صورت بازگشتی محاسبه می‌کند.
3. دو DFT به دست آمده را با استفاده از «ضرایب چرخش» (Twiddle factors) به صورت $e^{-2\pi i k/N}$ ترکیب می‌کند.

$$\begin{aligned} X[k] &= E[k] + e^{-2\pi i k/N} \, O[k] && \text{for } k = 0, \ldots, N/2 - 1 \\ X[k + N/2] &= E[k] - e^{-2\pi i k/N} \, O[k] && \text{for } k = 0, \ldots, N/2 - 1 \\[6pt] \text{where } E &= \text{DFT of even-indexed samples} \\ O &= \text{DFT of odd-indexed samples} \end{aligned}$$

این تقارن باعث می‌شود هر سطح از بازگشت به اندازه‌ی $O(N)$ کار ببرد و از آنجا که $\log_2(N)$ سطح داریم، کل محاسبات به $O(N \log N)$ می‌رسد.

الگوریتم FFT پایه‌دو (Radix-2) که در اینجا پیاده‌سازی شده، وقتی طول سیگنال توانی از 2 باشد بیشترین سرعت را دارد. هر بار که طول یک زیرمسئله فرد شود، الگوریتم به dft ساده پس‌روی (Fallback) می‌کند؛ به همین دلیل معمولاً سیگنال‌ها را با پدینگ صفر (Zero-padding) به توان بعدی 2 می‌رسانند تا محاسبات سریع‌تر شود. پیاده‌سازی‌های کتابخانه‌ای مانند numpy.fft با الگوریتم‌های Radix ترکیبی هر طولی را مدیریت می‌کنند، اما بیشترین بازدهی آن‌ها برای طول‌هایی است که فاکتورهای اول کوچک دارند.

تحلیل طیفی

طیف توان (Power spectrum) به صورت $\lvert X[k] \rvert^2$ تعریف می‌شود که مجذور اندازه‌ی هر ضریب فرکانسی است و نشان می‌دهد چه مقدار انرژی در هر فرکانس وجود دارد.

طیف فاز (Phase spectrum) برابر با $\angle X[k]$ است که آفست فاز هر فرکانس را مشخص می‌کند. در بیشتر کارهای تحلیلی، معمولاً طیف توان اهمیت دارد و از فاز صرف‌نظر می‌شود.

$$\begin{aligned} \text{Power at frequency } k: \quad P[k] &= |X[k]|^2 = \mathop{\mathrm{Re}}(X[k])^2 + \mathop{\mathrm{Im}}(X[k])^2 \\ \text{Phase at frequency } k: \quad \phi[k] &= \mathop{\mathrm{atan2}}(\mathop{\mathrm{Im}}(X[k]), \mathop{\mathrm{Re}}(X[k])) \end{aligned}$$

تفکیک‌پذیری فرکانسی

تفکیک‌پذیری فرکانسی در DFT به تعداد نمونه‌ها $N$ و نرخ نمونه‌برداری $f_s$ بستگی دارد.

$$\begin{aligned} \text{Frequency of bin } k: \quad & f_k = k \cdot f_s / N \\ \text{Frequency resolution:} \quad & \Delta f = f_s / N \\ \text{Maximum frequency:} \quad & f_{\max} = f_s / 2 \quad \text{(Nyquist)} \end{aligned}$$

برای تفکیک دو فرکانس نزدیک به هم، به نمونه‌های بیشتری نیاز دارید. برای ثبت فرکانس‌های بالا هم به نرخ نمونه‌برداری بالاتری نیاز خواهید داشت.

قضیه‌ی کانولوشن

این قضیه یکی از مهم‌ترین دستاوردهای پردازش سیگنال است و ارتباط مستقیمی با لایه‌های CNN دارد.

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

$$\begin{aligned} & x * h = \text{IFFT}(\text{FFT}(x) \cdot \text{FFT}(h)) \\[4pt] & \text{where } * \text{ is convolution and } \cdot \text{ is element-wise multiplication} \end{aligned}$$

دلیل اهمیت موضوع:

• کانولوشن مستقیم دو سیگنال با طول‌های $N$ و $M$ به $O(NM)$ عملیات نیاز دارد.
• کانولوشن مبتنی بر FFT تنها $O(N \log N)$ عملیات می‌برد: تبدیل هر دو، ضرب، و تبدیل معکوس به دامنه‌ی زمان.
• برای کرنل‌های بزرگ، کانولوشن با FFT به طرز چشمگیری سریع‌تر است.
• این سازوکار دقیقاً همان اتفاقی است که در لایه‌های کانولوشن با میدان دید (Receptive field) بزرگ رخ می‌دهد.

نکته: تبدیل DFT کانولوشن دایره‌ای (Circular convolution) را محاسبه می‌کند (سیگنال دور می‌زند و تکرار می‌شود). برای به دست آوردن کانولوشن خطی (بدون تکرار و دور زدن)، باید پیش از محاسبه، هر دو سیگنال را با اضافه کردن صفر به طول $N + M - 1$ برسانید.

پنجره‌گذاری

تبدیل DFT فرض می‌کند سیگنال متناوب است؛ یعنی تعداد $N$ نمونه را به عنوان یک دوره‌ی تناوب از سیگنالی بی‌نهایت تکرارشونده در نظر می‌گیرد. اگر سیگنال با یک مقدار شروع نشود و با همان به پایان نرسد، ناپیوستگی در مرزها به وجود می‌آید که خود را به شکل محتوای فرکانس‌بالای کاذب نشان می‌دهد. به این پدیده نشتی طیفی (Spectral leakage) می‌گویند.

پنجره‌گذاری (Windowing) با میل دادن تدریجی سیگنال به صفر در هر دو انتها پیش از محاسبه‌ی DFT، این نشتی را کم می‌کند.

پنجره‌های متداول:

$$\begin{aligned} \text{Hann window:} \quad & w[n] = 0.5 \left(1 - \cos\frac{2\pi n}{N-1}\right) \\ \text{Hamming window:} \quad & w[n] = 0.54 - 0.46 \cos\frac{2\pi n}{N-1} \end{aligned}$$

پنجره را با ضرب عنصربه‌عنصر در سیگنال پیش از محاسبه‌ی DFT اعمال کنید: X = DFT(x * w).

ویژگی‌های DFT

قضیه‌ی پارسوال بیان می‌کند که کل انرژی در هر دو دامنه یکسان است؛ یعنی انرژی در طول تبدیل پایسته می‌ماند.

ارتباط با کدگذاری‌های مکانی

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

$$\begin{aligned} PE(pos, 2i) &= \sin\left(pos / 10000^{2i/d_{\text{model}}}\right) \\ PE(pos, 2i+1) &= \cos\left(pos / 10000^{2i/d_{\text{model}}}\right) \end{aligned}$$

هر جفت بُعد $(2i, 2i+1)$ با فرکانس متفاوتی نوسان می‌کند. این فرکانس‌ها به صورت هندسی از فرکانس بالا (ابعاد 0 و 1) تا فرکانس پایین (ابعاد نهایی) چیده شده‌اند. این طراحی باعث می‌شود هر موقعیت، الگوی منحصر‌به‌فردی در تمام باندهای فرکانسی داشته باشد؛ شبیه به اینکه چطور ضرایب فوریه یک سیگنال را به صورت یکتا مشخص می‌کنند.

ویژگی‌های کلیدی حاصل از این ساختار:

• یکتایی: هیچ دو موقعیتی کدگذاری یکسانی ندارند.
• کران‌دار بودن مقادیر: توابع سینوس و کسینوس همیشه در بازه‌ی $[-1, 1]$ هستند.
• موقعیت نسبی: کدگذاری موقعیت $p+k$ را می‌توان به صورت تابعی خطی از کدگذاری موقعیت $p$ نمایش داد؛ در نتیجه مدل می‌تواند یاد بگیرد که به موقعیت‌های نسبی توجه (Attention) کند.

ارتباط با شبکه‌های عصبی کانولوشنی

یک لایه‌ی کانولوشن فیلتر (کرنل) آموزش‌دیده‌ای را با لغزاندن روی سیگنال یا تصویر اعمال می‌کند. از نظر ریاضی، این دقیقاً همان عمل کانولوشن است.

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

1. اعمال FFT روی ورودی
2. اعمال FFT روی کرنل
3. ضرب آن‌ها در دامنه‌ی فرکانس
4. اعمال IFFT روی نتیجه

پیاده‌سازی‌های استاندارد CNN از کانولوشن مستقیم استفاده می‌کنند (که برای کرنل‌های کوچک 3x3 سریع‌تر است). اما برای کرنل‌های بزرگ یا کانولوشن سراسری، روش‌های مبتنی بر FFT به مراتب سریع‌ترند. برخی معماری‌ها (مانند FNet) لایه‌ی اتنشن را کاملاً با FFT جایگزین کرده و به دقتی رقابتی با مرتبه‌ی زمانی $O(N \log N)$ به جای $O(N^2)$ دست یافته‌اند.

اسپکتروگرام‌ها و تبدیل فوریه‌ی کوتاه‌مدت

یک FFT منفرد، محتوای فرکانسی کل سیگنال را به شما می‌دهد، اما چیزی درباره‌ی زمان رخ دادن این فرکانس‌ها نمی‌گوید. یک سیگنال چیرپ (Chirp؛ سیگنالی که فرکانس آن در طول زمان افزایش می‌یابد) و یک آکورد (که تمام فرکانس‌هایش به طور هم‌زمان حضور دارند) می‌توانند طیف اندازه‌ی کاملاً یکسانی داشته باشند.

تبدیل فوریه‌ی کوتاه‌مدت (STFT) این مشکل را با محاسبه‌ی FFT روی پنجره‌های هم‌پوشان سیگنال حل می‌کند. نتیجه‌ی این کار یک اسپکتروگرام است: بازنمایی دوبعدی با زمان روی یک محور و فرکانس روی محور دیگر. شدت در هر نقطه نشان‌دهنده‌ی انرژی آن فرکانس در همان زمان است.

STFT procedure:
1. Choose a window size (e.g., 1024 samples)
2. Choose a hop size (e.g., 256 samples -- 75% overlap)
3. For each window position:
   a. Extract the windowed segment
   b. Apply a Hann/Hamming window
   c. Compute FFT
   d. Store the magnitude spectrum as one column of the spectrogram

اسپکتروگرام‌ها ورودی استاندارد مدل‌های یادگیری ماشین در حوزه‌ی پردازش صوت هستند. مدل‌های تشخیص گفتار (مانند Whisper یا DeepSpeech) روی مِل-اسپکتروگرام‌ها کار می‌کنند؛ یعنی اسپکتروگرام‌هایی که فرکانس‌های آن‌ها به مقیاس مل (Mel scale) نگاشت شده تا تطابق بهتری با ادراک شنوایی انسان از زیروبمی صدا داشته باشد.

الیاسینگ

اگر سیگنالی شامل فرکانس‌هایی بالاتر از $f_s/2$ (فرکانس نایکویست) باشد، نمونه‌برداری با نرخ $f_s$ نسخه‌های کاذبی از آن‌ها خواهد ساخت. یک سیگنال 90 هرتز که با نرخ 100 هرتز نمونه‌برداری شود، کاملاً شبیه به یک سیگنال 10 هرتز به نظر می‌رسد و هیچ راهی برای تمایز آن‌ها صرفاً از روی نمونه‌ها وجود ندارد.

Example:
  True signal: 90 Hz sine wave
  Sampling rate: 100 Hz
  Apparent frequency: 100 - 90 = 10 Hz

  The samples from the 90 Hz signal at 100 Hz sampling rate
  are identical to the samples from a 10 Hz signal.
  No amount of math can recover the original 90 Hz.

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

پدینگ صفر تفکیک‌پذیری را افزایش نمی‌دهد

یک تصور اشتباه رایج این است که اضافه کردن صفر (Zero-padding) پیش از اعمال FFT، تفکیک‌پذیری فرکانسی را بالا می‌برد؛ اما این‌طور نیست. پدینگ صفر تنها میان خانه‌های فرکانسی موجود درونیابی (Interpolation) می‌کند و طیف هموارتری تحویل می‌دهد، اما نمی‌تواند جزئیات فرکانسی جدیدی را که در نمونه‌های اصلی وجود نداشته آشکار کند.

تفکیک‌پذیری فرکانسی واقعی تنها به زمان مشاهده $T = N / f_s$ بستگی دارد. برای تفکیک دو فرکانس با فاصله‌ی $\Delta f$ از یکدیگر، حداقل به $T = 1 / \Delta f$ ثانیه داده نیاز دارید و هیچ میزانی از پدینگ صفر این محدودیت بنیادین را تغییر نمی‌دهد.

پیاده‌سازی از صفر

گام 1: پیاده‌سازی DFT از صفر

الگوریتم DFT با مرتبه‌ی $O(N^2)$ مستقیماً از روی تعریف پیاده‌سازی می‌شود.

import math

class Complex:
    ...

def dft(x):
    N = len(x)
    result = []
    for k in range(N):
        total = Complex(0, 0)
        for n in range(N):
            angle = -2 * math.pi * k * n / N
            w = Complex(math.cos(angle), math.sin(angle))
            xn = x[n] if isinstance(x[n], Complex) else Complex(x[n])
            total = total + xn * w
        result.append(total)
    return result

گام 2: تبدیل معکوس DFT

همان ساختار، با توان مثبت و تقسیم بر $N$.

def idft(X):
    N = len(X)
    result = []
    for n in range(N):
        total = Complex(0, 0)
        for k in range(N):
            angle = 2 * math.pi * k * n / N
            w = Complex(math.cos(angle), math.sin(angle))
            total = total + X[k] * w
        result.append(Complex(total.real / N, total.imag / N))
    return result

گام 3: پیاده‌سازی FFT (الگوریتم Cooley-Tukey)

الگوریتم بازگشتی FFT نیاز دارد که طول سیگنال حتماً توانی از 2 باشد. به دو بخش زوج و فرد تقسیم کنید، بازگشت بزنید، و با ضرایب چرخش ادغام نمایید.

def fft(x):
    N = len(x)
    if N <= 1:
        return [x[0] if isinstance(x[0], Complex) else Complex(x[0])]
    if N % 2 != 0:
        return dft(x)

    even = fft([x[i] for i in range(0, N, 2)])
    odd = fft([x[i] for i in range(1, N, 2)])

    result = [Complex(0)] * N
    for k in range(N // 2):
        angle = -2 * math.pi * k / N
        twiddle = Complex(math.cos(angle), math.sin(angle))
        t = twiddle * odd[k]
        result[k] = even[k] + t
        result[k + N // 2] = even[k] - t
    return result

گام 4: توابع کمکی تحلیل طیفی

def power_spectrum(X):
    return [xk.real ** 2 + xk.imag ** 2 for xk in X]

def convolve_fft(x, h):
    N = len(x) + len(h) - 1
    padded_N = 1
    while padded_N < N:
        padded_N *= 2

    x_padded = x + [0.0] * (padded_N - len(x))
    h_padded = h + [0.0] * (padded_N - len(h))

    X = fft(x_padded)
    H = fft(h_padded)

    Y = [xk * hk for xk, hk in zip(X, H)]

    y = idft(Y)
    return [y[n].real for n in range(N)]

استفاده از ابزارهای آماده

برای کارهای واقعی، از ماژول FFT کتابخانه‌ی NumPy استفاده کنید که متکی بر کتابخانه‌های C بسیار بهینه‌سازی‌شده است.

import numpy as np

signal = np.sin(2 * np.pi * 5 * np.arange(256) / 256)
spectrum = np.fft.fft(signal)
freqs = np.fft.fftfreq(256, d=1/256)

power = np.abs(spectrum) ** 2

positive_freqs = freqs[:len(freqs)//2]
positive_power = power[:len(power)//2]

برای پنجره‌گذاری و تحلیل‌های پیشرفته‌تر طیفی:

from scipy.signal import windows, stft

window = windows.hann(256)
windowed = signal * window
spectrum = np.fft.fft(windowed)

برای انجام کانولوشن:

from scipy.signal import fftconvolve

result = fftconvolve(signal, kernel, mode='full')

برای محاسبه‌ی اسپکتروگرام:

from scipy.signal import stft

frequencies, times, Zxx = stft(signal, fs=sample_rate, nperseg=256)
spectrogram = np.abs(Zxx) ** 2

ماتریس اسپکتروگرام ابعادی به صورت (n_frequencies, n_time_frames) دارد. هر ستون نشان‌دهنده‌ی طیف توان در یک پنجره‌ی زمانی است. مدل‌های یادگیری ماشین در حوزه‌ی صوت دقیقاً چنین داده‌ای را به عنوان ورودی مصرف می‌کنند.

خروجی نهایی

برای تولید فایل outputs/prompt-spectral-analyzer.md، اسکریپت code/fourier.py را اجرا کنید.

تمرین‌ها

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

2. اعتبارسنجی FFT در برابر DFT. سیگنالی تصادفی به طول 64 بسازید. هر دو تبدیل DFT ($O(N^2){}$) و FFT را محاسبه کنید. بررسی کنید که تمام ضرایب تا دقت 1e-10 هم‌خوانی داشته باشند. زمان اجرای هر دو تابع را روی سیگنال‌هایی به طول‌های 256، 512، 1024 و 2048 اندازه بگیرید. نسبت زمان اجرای DFT به FFT را ترسیم کنید.

3. اثبات تجربی قضیه‌ی کانولوشن. سیگنال x = [1, 2, 3, 4, 0, 0, 0, 0] و فیلتر h = [1, 1, 1, 0, 0, 0, 0, 0] را بسازید. کانولوشن دایره‌ای آن‌ها را یک بار مستقیم (با حلقه‌های تودرتو) محاسبه کنید؛ سپس آن را با FFT (تبدیل، ضرب، معکوس تبدیل) به دست آورید و مطابقت نتایج را بسنجید. حالا کانولوشن خطی را با پدینگ صفر مناسب پیاده کنید.

4. اثرات پنجره‌گذاری. سیگنالی بسازید که حاصل جمع دو موج سینوسی در فرکانس‌های 10 و 12 هرتز (بسیار نزدیک به هم) باشد. با نرخ 128 هرتز به مدت 1 ثانیه نمونه‌برداری کنید. طیف توان را بدون پنجره، با پنجره‌ی Hann و با پنجره‌ی Hamming محاسبه کنید. کدام پنجره تشخیص دو قله را راحت‌تر می‌کند؟ چرا؟

5. تحلیل کدگذاری مکانی. کدگذاری‌های مکانی سینوسی را برای d_model = 128 و max_pos = 512 تولید کنید. برای هر جفت موقعیت $(p_1, p_2)$، ضرب داخلی کدهای آن‌ها را محاسبه کنید. نشان دهید که حاصل‌ضرب داخلی فقط به $\lvert p_1 - p_2 \rvert$ بستگی دارد، نه به موقعیت‌های مطلق. با افزایش فاصله، چه بر سر ضرب داخلی می‌آید؟

راهنمایی‌ها

تمرین 1. با 128 نمونه در 1 ثانیه، خانه‌ی فرکانسی $k$ معادل $k$ هرتز است؛ بنابراین برای یک فرکانس صحیح، اندیس بزرگ‌ترین مقدار در خانه‌های 1 تا 63 همان فرکانس خواهد بود. یک عدد صحیح برای هرتز انتخاب کنید. اگر فرکانس بین دو خانه بیفتد، دامنه‌اش در خانه‌های مجاور پخش می‌شود و قله فقط نزدیک‌ترین خانه را نشان می‌دهد. از generate_signal و magnitude_spectrum استفاده کنید. پس از افزودن نویز، دقت کنید که قله همچنان در همان خانه بماند، هرچند سطح سایر خانه‌ها بالاتر می‌آید.

تمرین 2. طول سیگنال را توانی از 2 انتخاب کنید، چون fft سیگنال را به دو نیمه تقسیم می‌کند. مقدار max(abs(a - b)) را روی تمام ضرایب بررسی کنید. نسبت زمانی باید تقریباً با نرخ $N / \log_2 N$ رشد کند.

تمرین 3. برای کانولوشن دایره‌ای، مقدار $y[n] = \sum_k h[k] x[(n-k) \bmod 8]$ را به صورت دستی برای مقادیر $n$ برابر با 0، 1 و 2 محاسبه کرده و با سه مقدار نخستِ convolve_direct و convolve_fft مقایسه کنید. کانولوشن خطی طولی برابر با 6 = 1 - 3 + 4 دارد؛ پس ابتدا بررسی کنید که آیا طول 8 از قبل به اندازه‌ی کافی بزرگ هست یا خیر.

تمرین 4. از hann_window، hamming_window، apply_window و power_spectrum استفاده کنید و قله‌های تفکیک‌شده را بین 9 تا 13 هرتز بشمارید. چون 10 و 12 هرتز هر دو چرخه‌های کاملی را در 1 ثانیه طی می‌کنند، ممکن است حالت بدون پنجره بهترین نتیجه را بدهد. این کار را با 10.5 و 12.5 هرتز تکرار کنید تا نشتی طیفی‌ای را که پنجره‌ها کاهش می‌دهند مشاهده کنید.

تمرین 5. هر کدگذاری را با positional_encoding(pos, 128) محاسبه کنید. حاصل‌ضرب داخلی برای موقعیت‌های (0, 10) را با حاصل‌ضرب داخلی برای (100, 110) مقایسه نمایید؛ این دو مقدار باید تا دقت حدود 1e-10 یکسان باشند. سپس برای یک $p_1$ مشخص و ثابت، نمودار ضرب داخلی را بر حسب $\lvert p_1 - p_2 \rvert$ رسم کنید.

اصطلاحات کلیدی

برای مطالعه‌ی بیشتر

• Cooley & Tukey: An Algorithm for the Machine Calculation of Complex Fourier Series (1965) - مقاله‌ی بنیادین FFT که دنیای محاسبات را دگرگون کرد
• 3Blue1Brown: But what is the Fourier Transform? - بهترین معرفی بصری برای درک تبدیل‌های فوریه
• Lee-Thorp et al.: FNet: Mixing Tokens with Fourier Transforms (2021) - جایگزینی سلف‌اتنشن با FFT در ترنسفورمرها
• Smith: The Scientist and Engineer's Guide to Digital Signal Processing - کتاب درسی رایگان و جامعی که عمیقاً به مفاهیم FFT، پنجره‌گذاری و تحلیل طیفی می‌پردازد
• Vaswani et al.: Attention Is All You Need (2017) - کدگذاری‌های مکانی سینوسی برگرفته از تجزیه‌ی فرکانسی فوریه
• Radford et al.: Whisper (2022) - تشخیص گفتار با استفاده از مل-اسپکتروگرام‌ها به عنوان بازنمایی ورودی