前言

最近復現音樂驅動舞蹈的文章《Dancing-to-Music Character Animation》，用到了與傅里葉變換很相似的稱為常Q變換的方法去分割音樂，所以對傅里葉變換做了一個小了解，本文不深入各種亂糟糟的理論，比如什么蝶形算法啥的，單純討論離散傅里葉變換(DFT)，我們常見的是快速傅里葉變換(FFT)，其實FFT是對DFT的一個計算優化，主要是剔除DFT里面有周期性之類的被冗余計算，但是FFT的算法有點小復雜，有興趣深入理論的請移步如下幾篇博文：

如何理解和掌握快速傅里葉變換的計算和概念？

如何通俗地解釋什么是離散傅里葉變換？

傅里葉分析之掐死教程（完整版）更新于2014.06.06

傅里葉變換

快速傅里葉變換

第三章 離散傅里葉變換(DFT) 及其快速算法(FFT)

傅里葉級數和傅里葉變換是什么關系？

如何理解傅里葉變換公式？

An Introduction to the Fast Fourier Transform

FFT的詳細解釋，相信你看了就明白了

如果想仔細學習FFT，最好是看完上述推薦的博文并額外找資料，我是不想看了，因為我看著看著發現自己掉頭發了o(╯□╰)o

#簡介

傅里葉變換意思就是任何一個輸入信號都可以使用多個余弦波疊加而成，簡單的說就是把時序信號轉換成頻域信息。我們最終需要的就是找到這些余弦波的相關參數：幅值，相位。

復習一下三角函數的標準式：
$$y=Acos(wx+b)+k$$
$A$代表振幅，函數周期是$\frac{2\pi }{w}$，頻率是周期的倒數$\frac{w}{2\pi }$，$b$是函數初相位，$k$在信號處理中稱為直流分量。

在很多工具的實現中，余弦波的個數就是信號長度，但是在理論公式中會出現一個參數N，是采樣點，通常稱為N點FFT。
$$X(k)=\sum _{n=1}^{N}x(n)?\exp (?\sqrt{?1}?2\pi ?(k?1)?(n?1)/N),1<=k<=N.$$
上述公式就是DFT的求解方法了，不考慮它的優化方法FFT，我們直接在matlab中碼公式，最后與FFT的結果做對比即可驗證寫的算法對不對。

#實例

假設我們的輸入信號是函數是
$$S=0.2+0.7?\cos (2\pi ?50t+\frac{20}{180}\pi )+0.2?\cos (2\pi ?100t+\frac{70}{180}\pi )$$
可以發現直流分量是0.2，以及兩個余弦函數的疊加，余弦函數的幅值分別為0.7和0.2，頻率分別為50和100，初相位分別為20度和70度。

畫出信號圖：

Fs = 1000; % Sampling frequency T = 1/Fs; % Sampling period L = 1000; % Length of signal t = (0:L-1)*T; % Time vector S = 0.2+0.7*cos(2*pi*50*t+20/180*pi) + 0.2*cos(2*pi*100*t+70/180*pi) ; plot(1000*t(1:50),S(1:50)) title('Signal Corrupted with Zero-Mean Random Noise') xlabel('t (milliseconds)') ylabel('X(t)')

FFT變換

先使用matlab默認的快速傅里葉變換函數FFT求解疊加余弦的各參數

Y = fft(S); P2 = abs(Y/L); P1 = P2(1:L/2+1); P1(2:end-1) = 2*P1(2:end-1);

首先直接對原始信號進行傅里葉變換得到$Y$，它包含實部與虛部，然后求解歸一化后$Y$各項的模得到$P2$，由于matlab求解的結果包含對稱的兩個頻譜，稱為雙側頻譜，我們只需要取一半得到$P1$，此時只需要將除第一個元素外的元素乘以2即可得到幅值信息，其實求解的根本就是來源于$Y$，$Y$有多少項，就說明求解了多少個疊加的余弦波，接下來解釋如何求解各參數：

• 直流分量：就是Y的第一個值除以采樣頻率，一般來說是非復數
• 頻率：$\frac{采樣頻率}{采樣長度}?(第幾項?1)$，本例中采樣頻率是1000，長度也是1000，那么$Y$的第二項頻率就是1，第三項頻率是2，其實最終情況下，我們選取頻率不接近0的值。
• 幅值：$2?abs(\frac{Y各項}{采樣長度})$
• 初相位：$atan2(\frac{Y的虛部}{Y的實部})$轉角度制表示

從$P1$的圖中，我們很容易看出來幅值不接近0的位置分別是0，50，100附近，其實我們去看具體的位置發現是1,51,101，此三個位置的$Y$值分別為：

>> Y(1)ans =200.0000>> Y(51)ans =3.2889e+02 + 1.1971e+02i>> Y(101)ans =34.2020 +93.9693i

那么按照描述，我們得到：

• 直流分量：$\frac{Y(1)}{Fs}=0.2$

• 頻率：第51項和101項的頻率分別為50和100

• 幅值：$2?abs\left(\frac{Y(51)}{L}\right)=0.7$ 同理$2?abs\left(\frac{Y(101)}{L}\right)=0.2$

• 初相位：

  >> rad2deg(atan2(imag(Y(51)),real(Y(51))))ans =20.0000>> rad2deg(atan2(imag(Y(101)),real(Y(101))))ans =70.0000

【注1】: 在實際應用中，一般讓余弦波的的數量與信號長度相同，如果不相同，那就是理論中常說的N點FFT

【注2】：幅值是通過模求解的，但是模一般都是正數，如果疊加信號的幅值是復數呢？不要忘記了$?cos(x)=cos(x?\pi )$，也就是說如果幅值是負數，那么只需要在正數幅值的基礎上變化一初相位。比如把例子的函數換成：
$$S=0.2?0.7?\cos (2\pi ?50t+\frac{20}{180}\pi )+0.2?\cos (2\pi ?100t+\frac{70}{180}\pi )$$
那么求解頻率50對應余弦波的賦值為+0.7，初相位為-160

IFFT變換

顧名思義，IFFT就是逆傅里葉變換，用Y重構信號，其實我們通過Y都已經分析出了原始信號所需的余弦波的各參數，肯定能重構原始數據，這里還是做個實驗吧，用IFFT函數：

figure pred_X=ifft(Y); plot(pred_X,'r-') hold on plot(S,'b*')

DFT變換

按照公式手擼DFT，看看計算得到$Y$與matlab自帶的FFT結果是否一致
$$X(k)=\sum _{n=1}^{N}x(n)?\exp (?\sqrt{?1}?2\pi ?(k?1)?(n?1)/N),1<=k<=N.$$

%% DFT變換 % N % X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N. % n=1 DFT_X=zeros(1,L); N=L; for k=1:Lfor n=1:NDFT_X(k)=DFT_X(k)+S(n)*exp(-1j*2*pi*(k-1)*(n-1)/N);end end P2=abs(DFT_X/L); P1 = 2*P2(1:L/2+1); f = Fs*(0:(L/2))/L; figure plot(f,P1) xlabel('頻率') ylabel('振幅') title('DFT變換')

再計算與FFT求解的結果的誤差

%% FFT變換 FFT_X=fft(S); figure plot(abs(FFT_X-DFT_X)) title('DFT和FFT誤差')

IDFT

同樣按照公式手擼逆離散傅里葉變換
$$x(n)=\frac{1}{N}\sum _{k=1}^{N}X(k)?\exp (\sqrt{?1}?2\pi ?(k?1)?(n?1)/N),1<=n<=N.$$

%% IDFT變換 % N % x(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N. % k=1 DFT_rec_x=zeros(1,k); for n=1:Lfor k=1:LDFT_rec_x(n)=DFT_rec_x(n)+DFT_X(k)*exp( 1j*2*pi*(k-1)*(n-1)/N);end end DFT_rec_x=DFT_rec_x./N; rec_err=real(DFT_rec_x)-S; figure plot(rec_err) title('IFFT數據重構誤差')

與IFFT的結果對比一下：

%% IFFT變換 FFT_rec_x=ifft(FFT_X); figure plot(abs(DFT_rec_x-FFT_rec_x)) title('IDFT和IFFT誤差')

后記

折騰了這么多，其實就是為了一個字：懶。為了避免復雜的FFT的理論理解，我直接按照DFT的公式碼代碼，取得了與FFT一樣的結果，對于不喜歡復雜理論的同志，可以在代碼實現FFT的時候直接采用DFT的代碼作為替代品，雖然時間復雜度增大很多，但是理論理解起來簡單很多倍不是。

貼一下文章代碼，此外還有一個FFT的手動實現：鏈接：https://pan.baidu.com/s/1mUdclEQ3tgQUvZQ4XffN3g 密碼：08tg

等我不掉頭發了，再去糾結FFT的蝶形算法^_

博客已同步更新到微信公眾號中，有興趣可關注一波，微信公眾號與博客同步更新個人的學習筆記

總結

以上是生活随笔為你收集整理的【音频处理】离散傅里叶变换的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。