完整版教程下載地址：http://www.armbbs.cn/forum.php?mod=viewthread&tid=94547
第28章 FFT和IFFT的Matlab實現（幅頻響應和相頻響應）
本章主要講解fft，ifft和fftshift在matlab上的實現。
28.1 初學者重要提示
28.2 Matlab的FFT函數
28.3 Matlab的IFFT函數
28.4 Matlab的FFTSHIFT函數
28.5 總結
28.1 初學者重要提示
- 求解FFT相頻時的修正比較重要，本章做了一個簡易修正方法。重點看28.2.5小節。
28.2 Matlab的FFT函數
28.2.1 函數語法
Y = fft(x)
Y = fft(X,n)
Y = fft(X,n,dim)
28.2.2 函數定義
Y = fft(x) 和 y = ifft(X)分別用於實現正變換和逆變換，公式描述如下：
28.2.3 函數描述
Y = fft(X)
用快速傅里葉變換 (FFT) 算法計算 X 的離散傅里葉變換 (DFT)。
- 如果 X 是向量，則 fft(X) 返回該向量的傅里葉變換。
- 如果 X 是矩陣，則 fft(X) 將 X 的各列視為向量，並返回每列的傅里葉變換。
- 如果 X 是一個多維數組，則 fft(X) 將尺寸大小不等於 1 的第一個數組維度的值視為向量，並返回每個向量的傅里葉變換。
注意這里第一個尺寸不為1是指一個矩陣的第一個尺寸不為1的維。
比如一個矩陣是2*1，那么第一個尺寸不為1的維就是行（尺寸為2）。
X是 1*2*3表示第一個尺寸不為1的維就是列（尺寸為2）。
X為維數5*6*2的話，第一個尺寸不為1的維就是行（尺寸為5）。
Y = fft(X, n)
返回 n 點 DFT。如果未指定任何值，則 Y 的大小與 X 相同。
- 如果 X 是向量且 X 的長度小於 n，則為 X 補上尾零以達到長度 n。
- 如果 X 是向量且 X 的長度大於 n，則對 X 進行截斷以達到長度 n。
- 如果 X 是矩陣，則每列的處理與在向量情況下相同。
- 如果 X 為多維數組，則大小不等於 1 的第一個數組維度的處理與在向量情況下相同。
Y = fft(X, n, dim)
返回沿維度 dim 的傅里葉變換。例如，如果 X 是矩陣，則 fft(X,n,2) 返回每行的 n 點傅里葉變換。
28.2.4 FFT實例一：幅頻響應
傅里葉變換的一個常見用途就是查找埋藏在噪聲信號中的實際信號的頻率成分。下面我們考慮一個這樣的例子：
采樣率是1000Hz ，信號由如下三個波形組成。
（1）50Hz的正弦波、振幅0,7。
（2）70Hz正弦波、振幅1。
（3）均值為0的隨機噪聲。
實際運行代碼如下：
Fs = 1000; %采樣率 T = 1/Fs; %采樣時間單位 L = 1000; %信號長度 t = (0:L-1)*T; %時間序列 x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); %原始信號 y = x + 2*randn(size(t)); %原始信號疊加了噪聲后 plot(Fs*t(1:50),y(1:50)); %繪制波形 title('原始信號+零均值隨機噪聲 '); xlabel('時間單位:ms');
運行Matlab后，顯示波形如下：
通過上面的截圖，我們是很難發現波形中的頻率成分，下面我們通過FFT變換，從頻域觀察就很方便了，Matlab運行代碼如下：
Fs = 1000; %采樣率 T = 1/Fs; %采樣時間單位 L = 1000; %信號長度 t = (0:L-1)*T; %時間序列 x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); %原始信號 y = x + 2*randn(size(t)); %原始信號疊加了噪聲后 NFFT = 2^nextpow2(L); %求得最接近采樣點的2^n，由於上面是1000點，那么最近的就是1024點。 Y=fft(y,NFFT)/L; % 進行FFT變換，除以總的采樣點數，方便觀察實際值。 f = Fs/2*linspace(0,1,NFFT/2+1); %頻率軸，這里只顯示Fs/2部分，另一半是對稱的。 plot(f,2*abs(Y(1:NFFT/2+1))) %繪制波形 title('幅頻相應'); xlabel('頻率'); ylabel('幅度');
從上面的幅頻相應，我們可以看出，兩個正弦波的頻譜並不是准確的0.7和1，而是比較接近，這個就是我們在上節教程中所示的頻譜泄露以及噪聲的干擾。
28.2.5 FFT實例二：相頻響應（重要）
這里我們以采樣兩個余弦波組成的信號為例進行說明，並求出其幅頻和相頻響應。
（1）50Hz的余弦波，初始相位60°，振幅1.5。
（2）90Hz的余弦波、初始相位60°，振幅1。
（3）采樣率256Hz，采集256個點。
Matlab上面運行的代碼如下：
Fs = 256; % 采樣率 N = 256; % 采樣點數 n = 0:N-1; % 采樣序列 t = 0:1/Fs:1-1/Fs; % 時間序列 f = n * Fs / N; %真實的頻率 x = 1.5*cos(2*pi*50*t+pi/3)+ cos(2*pi*90*t+pi/3) ; %原始信號 y = fft(x, N); %對原始信號做FFT變換 Mag = abs(y); %求FFT轉換結果的模值 subplot(2,1,1); plot(f, Mag); %繪制幅頻相應曲線 title('幅頻相應'); xlabel('頻率/Hz'); ylabel('幅度'); subplot(2,1,2); plot(f, angle(y)*180/pi); %繪制相頻響應曲線，注意這將弧度轉換成了角度 title('相頻響應方式一'); xlabel('頻率/Hz'); ylabel('相角');
運行后求出的幅頻相應和相頻響應結果如下：
求出的幅頻響應沒問題，而相頻響應雜亂無章，造成這個問題的根本原因很多頻段的幅值非常小，他們的相角可以不顯示出來，這樣就可以方便的查看相頻響應了。基於此，修正后的代碼如下：
Fs = 256; % 采樣率 N = 256; % 采樣點數 n = 0:N-1; % 采樣序列 t = 0:1/Fs:1-1/Fs; % 時間序列 f = n * Fs / N; %真實的頻率 x = 1.5*cos(2*pi*50*t+pi/3)+ cos(2*pi*90*t+pi/3) ; %原始信號 y = fft(x, N); %對原始信號做FFT變換 Mag = abs(y); %求FFT轉換結果的模值 subplot(3,1,1); plot(f, Mag); %繪制幅頻相應曲線 title('幅頻相應'); xlabel('頻率/Hz'); ylabel('幅度'); subplot(3,1,2); plot(f, angle(y)*180/pi); %繪制相頻響應曲線，注意這將弧度轉換成了角度 title('相頻響應方式一'); xlabel('頻率/Hz'); ylabel('相角'); subplot(3,1,3); plot(f, angle(y)*180/pi.*(Mag>=100)); %繪制相頻響應曲線，注意這將弧度轉換成了角度 title('相頻響應方式二'); xlabel('頻率/Hz'); ylabel('相角');
上面代碼中Mag>=100是關鍵，僅展示FFT后幅值大於100的相角。下面再來看Matlab的效果：
可以看到已經完全沒問題了，求出了頻率50Hz的余弦初相為60°左右，頻率90Hz的余弦初相也是60°。
28.3 Matlab的IFFT函數
28.3.1 函數語法
y = ifft(X)
y = ifft(X,n)
y = ifft(X,[],dim)
y = ifft(X,n,dim)
y = ifft(..., 'symmetric')
y = ifft(..., 'nonsymmetric')
28.3.2 函數描述
y = ifft(X)
此函數用於返回向量X的離散傅立葉變換（DFT）逆變換結果，計算時使用快速傅里葉算法（Fast Fourier transform (FFT)）。
y = ifft(X,n)
此函數用於返回n點的IDFT。
y = ifft(X,[],dim)
y = ifft(X,n,dim)
上面兩個函數用於實現指定維度的IFFT運算。
28.3.3 IFFT實例
下面我們對信號：0.7*sin(2*pi*50*t) + sin(2*pi*120*t)求FFT和IFFT，並繪制原始信號和轉換后的信號。Matlab上運行的代碼如下：
Fs = 1000; %采樣率 T = 1/Fs; % 采樣時間 L = 1024; % 信號長度 t = (0:L-1)*T; % 時間序列 y = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); %50Hz正弦波和120Hz的正弦波的疊加 subplot(2,1,1); plot(Fs*t(1:50), y(1:50)); %繪制原始信號 title('原始信號'); Y = fft(y, L); %分別調用正變換和逆變換 Z = ifft(Y); subplot(2,1,2); plot(Fs*t(1:50), Z(1:50)); %繪制逆變換后的波形 title('FFT和IFFT轉換后的信號');
運行后求出的結果如下：
通過上面的運行結果可以看出，轉換后的波形與原始的波形基本是一樣的。
28.4 Matlab的FFTSHIFT函數
fftshift的作用正是讓正半軸部分和負半軸部分的圖像分別關於各自的中心對稱。因為直接用fft得出的數據與頻率不是對應的，fftshift可以糾正過來
以下是Matlab的幫助文件中對fftshift的說明：
Y = fftshift(X) rearranges the outputs of fft, fft2, and fftn by moving the zero-frequency component to the center of the array. It is useful for visualizing a Fourier transform with the zero-frequency component in the middle of the spectrum. For vectors, fftshift(X) swaps the left and right halves of X。
下面我們在Matlab上面實現一個如下的代碼來說明fftshift的使用：
Fs = 256; % 采樣率 N = 256; % 采樣點數 n = 0:N-1; % 采樣序列 t = 0:1/Fs:1-1/Fs; % 時間序列 f = (-N/2:N/2-1) * Fs / N; %真實的頻率 x = 1.5*sin(2*pi*50*t+pi/3) ; %原始信號 y = fft(x, N); %對原始信號做FFT變換 Mag = abs(y); %求FFT轉換結果的模值 subplot(2,1,1); plot(f, Mag); %繪制幅頻相應曲線 title('fft幅頻相應'); xlabel('頻率/Hz'); ylabel('幅度'); z = fftshift(y); %對FFT轉換后的結果做偏移。 subplot(2,1,2); plot(f, z); %繪制幅頻相應曲線 title('fftshift幅頻相應'); xlabel('頻率/Hz'); ylabel('幅度');
Matlab的運行結果如下：
通過上面的運行結果我們可以看到，經過fftshift的調節后，正弦波的中心頻率正好對應在了相應的50Hz頻率點。使用fftshift還有很多其它的好處，有興趣的可以查找相關的資料進行了解。
28.5 總結
本章節主要講解了fft，iff和fftshift的基本用法，如果要深入了解，一定要多練習，多查資料和翻閱相關書籍。