『壹』 音頻演算法入門-傅里葉變換
上一篇文章中講了一個時域處理的演算法wsola，接下來會學習頻域處理演算法，在這之前必須得對頻域有所了解，這就不得不提傅里葉變換了，本文的目的是讓大家學會用傅里葉變換公式和傅里葉逆變換公式進行計算。數學公式是人們對世界中的現象的描述，我們學習數學公式也不該只停留在使用公式來解決問題的層次，得明白公式到底在描述什麼現象，從這些天才數學家的角度來看世界。懂的地方可跳過。項目地址在文章末尾給出。
我直接說結論，傅里葉級數公式包含了傅里葉變換和傅里葉逆變換（不嚴謹的說就是這么回事）。
先簡單說下具體關系，法國數學家傅里葉發現，任何周期函數都可以用正弦函數和餘弦函數構成的無窮級數來表示，這種表示方式就是傅里葉級數。假如有個波形比較復雜的周期函數，那麼找出能用來構成這個周期函數的正弦函數和餘弦函數的頻率的方法就叫做傅里葉變換，用這些頻率的正弦函數和餘弦函數疊加起來表示這個周期函數的方法就叫做傅里葉逆變換。
再從公式中看下他們的關系，首先介紹傅里葉級數到底是什麼，首先級數是指將數列的項依次用加號連接起來的函數。這么說可能大家還不理解，舉個例子：e^x=1+x/1!+x^2/2!+...x^n/n!....，等號左邊是指數函數，等號右邊就是級數。傅里葉級數公式如下：
我們主要看這個指數形式的傅里葉級數公式，把求和符號去掉，展開一下就是f(t)=Fa*e^jaω0t+Fb*e^jbω0t+Fc*e^jcω0t+Fd*e^jdω0.....。現在看下面的周期函數疊加效果圖，圖中顯示的是3個周期函數分別在坐標軸（橫軸時間，縱軸幅度）的圖像，寫成傅里葉級數形式就是f(t)=fa(t)+fb(t)+0+0....，這就是傅里葉級數公式要描述的現象。其中Fa*e^jaω0t=fa(t),Fb*e^jbω0t=fb(t),Fc*e^jcω0t=0....。
看下圖的傅里葉變換和逆變換公式，你會發現傅里葉逆變換公式和傅里葉級數公式極其相似，而傅里葉級數系數公式Fn又和傅里葉變換公式極其相似。所以對一個周期函數進行傅里葉級數展開的過程可以認為是先做傅里葉變換再做傅里葉逆變換的過程。
上圖就是傅里葉變換公式也叫連續傅里葉變換公式，有個很重要的事情，就是傅里葉變換公式和逆變換公式一定要一起給出，不然就會讓人誤解，你們在網上會看到各種各樣的寫法，但這些寫法都是對的，常見的如下圖所示。
為了方便後面的講解我把角頻率ω換成2πf，如上圖所示，ω是希臘字母讀作Omega，大寫是Ω，小寫是ω，以後這兩個字母會經常看到，都是等於2πf。不要和電學中的電阻單位搞混了，要明白字母只不過是一個符號而已，在不同學科領域都是混著用的，只要不和自己公式中其他字母沖突就行，例如上圖傅里葉變換公式中的j其實就是虛數單位i，一般時候我們會把虛數單位寫成i，但因為傅立葉變換經常用於電學解決一些問題，為了不和電流符號i混淆，所以公式就把i寫成j 。
要想了解傅里葉變換公式，首先要了解歐拉公式e^ix=cosx+isinx在圖像中的含義。以實部的值cosx作為橫坐標值，虛部sinx的值作為縱坐標值，x的取值從負無窮到正無窮，畫出所有的e^ix點後，你會發現這些點會形成一個周期為2π的圓。如下圖1所示（如果不理解，建議看3Blue1Brown的視頻，視頻連接：https://www.bilibili.com/video/BV1pW411J7s8）
所以歐拉公式e^ix其實就是隨著x的增大而在坐標繫上逆時針畫圓的過程，那麼e^-ix就表示順時針畫圓，e^-i2πx就表示畫圓的速度提高2π倍，也就是說x從0到1的過程就是順時針畫出一個完整圓的過程（當然x從1到2或者2到3等等，都能畫出一個完整的圓），把x換成t後，e^-i2πt表示每秒都會順時針畫出一個圓。e^-i2πft表示每秒都會順時針畫出f個圓。f(t)表示t時刻的振幅，f(t)函數畫出來就是時域波形圖。f(t)*e^-i2πft表示每經過1秒會順時針畫出f個圓，並在畫圓的同時，t時刻的圓半徑要乘上t時刻的振幅，其實就是以每秒的音頻振幅數據繞f圈的速度進行旋轉纏繞（為了方便理解，沒有用復雜的音頻數據，用的是一個頻率為3的正弦波音頻做的實驗，請看下圖2，圖的上半部分是時域波形圖，圖的左下角是f等於0.4的時候，用公式f(t)*e^-i2πft在實部和虛部構成的坐標系畫的圖，圖的右下角是頻譜圖，頻譜圖的橫坐標是頻率，縱坐標是振幅，振幅的值就是左下角圖中數據形成的圖案的質心（圖中的紅點）到坐標系原點的距離的2倍）。當改變f的值，你會發現數據大多數時候是和我們想的一樣，以坐標系原點為圓心環繞著，也就是振幅一直都是0，但是當f的值，也就每秒的圈數等於該音頻數據的頻率時，你會發現一個神奇的現象，那就是所有的數據會在實部或虛部坐標軸的一側形成一個圓（如下圖3所示，如此一來就知道這段音頻數據包含了一個頻率為3振幅為0.5的正弦波）。所以將多個正弦波疊加的音頻數據用傅里葉公式，f從負無窮到正無窮遍歷一遍，就可以把這個音頻數據里包含的正弦波都一一找出來。（如果不理解，建議看3Blue1Brown的視頻，視頻連接：https://www.bilibili.com/video/BV1pW411J7s8）
平時我們說的對音頻進行傅里葉變換處理，其實說的是短時離散傅里葉變換。短時離散傅里葉變換的公式（也可以直接叫做離散傅里葉變換公式）如下。
下面將教大家如何理解這個公式。上面說的連續傅里葉變換公式中有兩個原因導致我們無法使用，第一點要求是音頻數據的時間從負無窮到正無窮，第二點要求是任意時間t都要有幅度值x(t)才能代入公式進行計算。所以為了解決這兩個問題，把公式變為短時且離散的傅里葉變換公式，這個公式可以把一段時間（時間假設為Ts秒）的離散音頻數據（有N個采樣數據）進行傅里葉變換。你可以把離散傅里葉變換公式理解成連續傅里葉變換的變形，最重要的一點是連續傅里葉變換公式的f和離散傅里葉變換公式的k不是一個意思，他們的關系是k=f*Ts。所以離散傅里葉變換公式也可以寫成F(f)=1/n*∑f(t)*e^-j2πf*Ts*n/N，其中的Ts*n/N對應的就是連續傅里葉變換公式的t，只不過這個t沒辦法取任意時間了，t的取值也就隨著n的取值成為了離散的時間點，所以前面的系數由1/2π變為1/N。這樣這兩個公式就對應起來了。下面將進一步詳細介紹這個公式。
上一段說了k=f*Ts，這段我來解釋下為什麼，其實離散傅里葉變換公式中k表示的是這段Ts秒的音頻數據環繞坐標系原點的圈數，所以k並不是連續傅里葉變換公式里的頻率f，而頻率f指的是1秒鍾震盪的次數，在這個公式中頻率f也對應著1秒的音頻數據環繞的圈數，所以真正的頻率f=k/Ts。
有人可能會好奇，那為什麼不把離散傅里葉變換公式的自變數k換成f呢，這樣不是更好理解嗎？是會更好理解，但是沒有必要，用f的話還要做一次無用的換算。因為采樣點只有N個的原因，k的取值范圍就被限制住了，k的取值范圍只能是0~N-1的整數，這也是為什麼用k來做自變數而不是用f的原因。
還有人可能會好奇，傅里葉逆變換到底是怎麼把頻域的信息還原回時域的，其實公式計算出來的F(k)是一個復數，這個復數包含了這個頻率的周期函數的振幅和相位的信息，假設F(k)=a+ib，，F(k)的模|F(k)|=(a^2+b^2)^1/2，頻率f=k/Ts時的振幅為|F(k)|*2（因為求出來的值相當於圓心，但實際上振幅是圓離圓心最遠點到坐標原點的距離，所以要乘2），頻率f=k/Ts時的相位為arctan(b/a)。所以如果你知道一個周期函數包含了哪些頻率的周期函數，並且你這到這些周期函數的振幅和相位，你就可以像下圖一樣把fa(t)和fb(t)疊加在一起還原回f(t)。傅里葉逆變換的做法略有不同，但意思就是這么個意思，理解了離散傅里葉變換公式的計算，逆變換其實也是差不多代入數值計算就是了。（如果不理解怎麼用離散傅里葉變換公式計算，建議看視頻，視頻里有離散傅里葉變換完整的計算過程，視頻連接：https://www.hu.com/zvideo/1276595628009377792）
快速傅里葉變換推薦看下面兩個視頻
https://www.bilibili.com/video/BV1za411F76U
https://www.bilibili.com/video/BV1Jh411d7CN
下面是我用java實現的離散傅里葉變換及逆變換和快速傅里葉變換及逆變換，從他們的運行時間就可以看出來快速傅里葉變換快得多。（學完快速傅里葉變換再想想頻譜為何Y軸對稱？為何N/2對稱？）
『貳』 如何使用fft函數進行編程序和進行快速傅里葉逆變換
在圖象處理的廣泛應用領域中，傅立葉變換起著非常重要的作用，具體表現在包括圖象分析、圖象增強及圖象壓縮等方面。
fftshift的作用正是讓正半軸部分和負半軸部分的圖像分別關於各自的中心對稱。因為直接用fft得出的數據與頻率不是對應的，fftshift可以糾正過來。
假設f（x,y）是一個離散空間中的二維函數，則該函數的二維傅立葉變換的定義如下：
p=0,1…M-1 q=0,1…N-1 (1)
或 p=0,1…M-1 q=0,1…N-1 （2）
離散傅立葉反變換的定義如下：
m=0,1…M-1 n=0,1…N-1（3）
F（p,q）稱為f（m,n）的離散傅立葉變換系數。這個式子表明，函數f（m,n）可以用無數個不同頻率的復指數信號和表示，而在頻率（w1，w2）處的復指數信號的幅度和相位是F（w1，w2）。
2、MATLAB提供的快速傅立葉變換函數
（1）fft2
fft2函數用於計算二維快速傅立葉變換，其語法格式為：
B = fft2(I)
B = fft2(I)返回圖象I的二維fft變換矩陣，輸入圖象I和輸出圖象B大小相同。
例如，計算圖象的二維傅立葉變換，並顯示其幅值的結果，如圖所示，其命令格式如下
load imdemos saturn2
imshow(saturn2)
B = fftshift(fft2(saturn2));
imshow(log(abs(B)),[],'notruesize')
（2）fftshift
MATLAB提供的fftshift函數用於將變換後的圖象頻譜中心從矩陣的原點移到矩陣的中心，其語法格式為：
B = fftshift(I)
對於矩陣I，B = fftshift(I)將I的一、三象限和二、四象限進行互換。
（2）ifft2
ifft2函數用於計算圖象的二維傅立葉反變換，其語法格式為：
B = ifft2(I)
B = ifft2(A)返回圖象I的二維傅立葉反變換矩陣，輸入圖象I和輸出圖象B大小相同。其語法格式含義與fft2函數的語法格式相同，可以參考fft2函數的說明。
如果信號是二維的，用上面的函數即可！直接調用。
如果信號是一維的，給你下面的例子，你應該能明白！
clear
fs=100;N=128; %采樣頻率和數據點數
n=0:N-1;t=n/fs; %時間序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信號
y=fft(x,N); %對信號進行快速Fourier變換，逆變換函數為ifft
mag=abs(y); %求得Fourier變換後的振幅
f=n*fs/N; %頻率序列
subplot(2,2,1),plot(f,mag); %繪出隨頻率變化的振幅
xlabel('頻率/Hz');
ylabel('振幅');title('N=128');grid on;
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %繪出Nyquist頻率之前隨頻率變化的振幅
xlabel('頻率/Hz');
ylabel('振幅');title('N=128');grid on;
%對信號采樣數據為1024點的處理
fs=100;N=1024;n=0:N-1;t=n/fs;
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信號
y=fft(x,N); %對信號進行快速Fourier變換
mag=abs(y); %求取Fourier變換的振幅
f=n*fs/N;
subplot(2,2,3),plot(f,mag); %繪出隨頻率變化的振幅
xlabel('頻率/Hz');
ylabel('振幅');title('N=1024');grid on;
subplot(2,2,4)
plot(f(1:N/2),mag(1:N/2)); %繪出Nyquist頻率之前隨頻率變化的振幅
xlabel('頻率/Hz');
ylabel('振幅');title('N=1024');grid on;