學號：17020120039 姓名：韋濤 轉載自：https://blog.csdn.net/qq_37335890/article/details/85040277? 請求和清風明月? 《離散傅里葉變換DFT與FFT，MATLAB的FFT函數(shù)使用（原創(chuàng)）——如何使用fft（）繪制出真正的頻譜圖像》
以前一直對MATLAB中fft（）函數(shù)的使用一直存在疑惑，為什么要加一
些參數(shù)，并且如何確定這些參數(shù)，也查了許多資料，但很多都感覺只是
表面一說根本沒有講清其本質(zhì)。但隨著學習的推進，慢慢有所領悟，所
以打算把自己的一些所懂分享下，有什么問題也希望大家指正。
本文主要先對DFT、FFT的一些概念進行介紹，然后通過MATLAB仿真
進行fft（）分析，從而解釋上述參數(shù)
一、DFT與FFT
首先是對DFT與FFT的一些概念上的介紹，其實FFT與DFT是等價的，他們實現(xiàn)的功能是一樣的，只是FFT是DFT的算法優(yōu)化，因為畢竟要用電腦來計算，DFT算的太慢了，就優(yōu)化下也就成了FFT。所以此處我們對DFT與FFT的介紹是等價的。
那么我們就來介紹DFT，它也被叫做離散傅里葉變換，其實它就是DFS離散傅里葉級數(shù)的時域頻域主值序列，或者也是DTFT離散時間傅里葉變換的頻域采樣。至于DFS與DTFT相信大家也是明白，那么很多人好奇為什么還要DFT這玩意，主要呢，還是因為計算機，因為計算機不可能處理無限長的信號而DFS和DTFT要么時域無限長，要么頻域無限長，所以就搞了個DFT。
好的，那么我們就來直接看公式（這里我就假想大家都對其來源及原因都清楚了）：
注意看上圖，這就是我們接下來的所有依據(jù)，X（k）是正變換，x（n）反變換，此處務必注意x（n）那個式子有個1/N，這對后面的理解很關鍵。
好的那么接下來我們就可以利用MATLAB來進行分析了。
二、MATLAB的FFT函數(shù)使用
首先說明下什么情況下我們要用FFT，這是很簡單的，但還是要說說：因為現(xiàn)實世界都是連續(xù)的信號，相信我們要分析它時單純靠我們?nèi)怂闶呛苈闊┑?，所有這些都是需要計算機進行計算，那么答案就有了，我們其實都是處理連續(xù)信號，也就是說我們都是想要FFT來分析連續(xù)信號。
OK
1、那么這樣FFT第一個步驟就出現(xiàn)了，那就是A/D轉換，也就是對連續(xù)信號采樣。而這里我們就要確定一個參數(shù)采樣頻率Fs，說到采樣，大家肯定會想到采樣定理，其內(nèi)容就是當對信號進行以采樣頻率為Fs>=2fc的采樣時，信息不急丟失（fc為原始連續(xù)信號的最大截止頻率），那么這樣第一個參數(shù)Fs就這樣確定了，其應該滿足Fs>=2fc，這時采樣周期為Ts=1/Fs。
2、既然確定了采樣頻率Fs，那么我們就要想信號的時域應該怎么確定，換句話說，我們將采樣好的信號給計算機處理時，應該給計算機多少個。這里就出現(xiàn)了第二個參數(shù)即序列個數(shù)L，好的，在這里由于第一個參數(shù)確定的Fs，相應的也是采樣周期Ts，所以我們就可以得出時域信號的橫坐標時間的變化范圍，也就是
t=（0：L-1）Ts
那么這里就確定了N嗎？不，我們這里還沒確定，對于L的確立我們需要上面這個語句的幫助。那么究竟怎么確定呢，這里我們就要想到DFT的來源了：為什么DFT可以用來進行表示信號的頻譜響應，因為其時域頻域都是信號的相應時域頻域的主值區(qū)間，所以說我們要用DTF可以表示這個信號時，應盡可能的包含這個信號的一個周期或整數(shù)倍周期也就是上述t的范圍應該是信號的一個周期范圍或整數(shù)倍周期范圍，否則則會發(fā)生頻譜泄露，至于頻譜泄露其實就是出現(xiàn)了本來沒有的頻率分量。比如說，10Hz的cos（2π10t），本來只有一種頻率分量f=10Hz，但分析結果卻包含了與10Hz頻率相近的其它頻率分量。下述將會進行仿真分析。所以說此處L的確定應是LxTs=mT（T為信號的最小周期，m為正整數(shù)），常常m=1。那么或許有人問怎么平常我看到的都不是這樣，L都是取和Fs一樣的值（注意為了頻率分辨率為1Hz，下文會討論）。其實對于這個也有它的道理，不妨我們推到下：
LxTs=L/Fs=1，也就是說其表達的時域t范圍為1，那么平常我們的信號最小周期一般都小于1的，要是當信號最小周期大于1時其就是錯誤的啦，所以這個只是對大部分來說是有用的。
那么說了這么多究竟怎么確定，要是按LxTs=mT確定會不會太麻煩，確實有時我們真的不能知道確定的信號最小周期，所以這里就放寬點，只要保證LxTs>=T
就可以，當然如果可以取整數(shù)倍就最好了。
OK
可是上面說又會發(fā)生頻率泄露，這可怎么辦？那么這里就只能盡可能提高N-DFT的N點啦，這樣能使盡可能更多的能量落到正確的頻率點上，也就是說這樣誤差會更小，這樣就能達到我們的目的。
3、于是第三個參數(shù)N-DFT的N就出現(xiàn)了，根據(jù)上述2情況，我們知道N的提高有利于提高精確的，那么它還有其它約束嗎？答案肯定有的，這里就需要些理論知識了：
先來了解下頻域采樣定理：頻域以N點為采樣周期的采樣，在時域就是原序列的以周期N的延拓。
那么這樣也就是說上述2中的序列長度L必須小于或等于N，否則就會發(fā)生時域混疊，所以說L>=N這就是限制條件，常常我們?nèi)=N，當發(fā)生頻譜泄露比較嚴重時我們就可以考慮將L增大，這樣就可以減小頻譜泄露的的影響。
綜上，以上就是MATLAB相關的所有參數(shù)了，它們分別就采樣頻率Fs，原序列長度L，N-DFT變換的N，其主要確定關系總結如下：
Fs>=2fc（fc為信號最高截止頻率）
LxTs>=T（T為信號的最小周期，若能取LxTs=T則應盡可能取，可以避免頻譜泄露）
N>=L（N常常應為2^m，因為FFT運算就是不斷模2進行DFT，所以N為2的次方利于提高速度）
對于其常取值，F(xiàn)s=1024，N=1024，L=1024，注意滿足上述條件就好
說完了這些參數(shù)，我們就講講MATLAB的fft函數(shù)，其用法主要Y=fft（x，N），x為原信號序列，Y為DFT（也就是FFT）變換后的，但是大家會發(fā)現(xiàn)其變換后幅度變了而且變了很大，其實就是跟最初的那個式子有關，其頻域幅值增大了N倍。
并且發(fā)現(xiàn)其橫軸也都是整數(shù)，所以這時也要對其橫軸坐標進行變換，即要乘以相應的頻率分辨率Fs/N，所以說上面為什么，常常也取Fs=L，N=L,這樣為了頻率分辨率為1Hz，對于這些我們需用以下MATLAB語句進行實現(xiàn)
Y=fft(x,N);%N一定到大于信號序列x的長度，不過一般等于，決不能小于
%因為若要小于時域就會發(fā)生混疊
Y=abs(Y);
Y=Y./N;
%Y=Y.*2./N;Y(1)=Y(1)./2;若考慮將負頻率的幅度折算到正頻率時應這樣處理
%因為變換后有個N的乘積因子的影響，根據(jù)DFT公式可知，故消除其影響
f=(0:N-1)*Fs./N;
%頻率正?；驗樽儞Q后橫坐標是每個點對應的個數(shù)，故應轉化成實際的頻率
%由于頻譜是對稱的，且周期的故常常只畫一半如下
subplot(2,1,1);
plot(f(1:N/2),Y(1:N./2));
%當然也可以畫平常我們見到的有對稱的如下
f=f-Fs./2;
subplot(212);
plot(f,fftshift(Y));
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
上述代碼也有相關解釋，那么我也就對相關幾個語句進行介紹解釋：
f=(0:N-1)*Fs./N;這是N-DFT變換后的橫軸坐標，也就是真實頻率值；
Y=Y./N;這是對DFT變換后幅值處理，也就是要除于N才是真實的幅值
%Y=Y.*2./N;Y(1)=Y(1)./2;而這個是鑒于無真實的負頻率，其與正頻率對稱的原理，把負頻率去掉，全變?yōu)檎l率的幅值處理。
下面語句是為了變成我們平?？吹接姓擃l率的圖像
f=f-Fs./2;
subplot(212);
plot(f,fftshift(Y));
三、仿真試驗
OK接下來就是仿真試驗啦。
我們直接考慮x（t）=cos（2π10t）；其周期T=0.1s，fc=10Hz
1、首先考慮Fs=512，L=512，N=512，滿足上述三個關系吧，
Fs=512>=2fc；LxTs=1s>=T（也是10T，故沒頻率泄露）；N=512>=L=512；那么我們看下仿真結果是否對得到：
可以看到只有一個頻率，可以完全表示，故與上述推理一樣。代碼如下
Fs=512;%采樣頻率
Ts=1/Fs;%采樣周期
N=512;%N—DFT
L=512;%原信號序列長度
t=(0:L-1).*Ts;%時域自變量
x=cos(10*2*pi*t);%原信號
subplot(311);
plot(t,x);
title("原信號")
xlabel("t/s")
grid on
Y=fft(x,N);%fft變換
Y=abs(Y)./N;%實際幅值變換
f=(0:N-1)*Fs./N;%實際頻率變換
subplot(312);
plot(f(1:N/2),Y(1:N./2));
title("N-DFT變換幅頻響應單邊")
xlabel("f/Hz")
grid on
f=f-Fs./2;%移位
subplot(313);
plot(f,fftshift(Y));%移位
title("N-DFT變換幅頻響應雙邊")
xlabel("f/Hz")
grid on
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
2、那么我們再來考慮Fs=512，L=128，N=128，滿足上述三個關系吧，
Fs=512>=2fc；LxTs=0.15s>=T（不是T的整數(shù)倍，故有頻率泄露）；N=128>=L=128；那么我們看下仿真結果是否對得到（這次頻譜用stem繪制更明顯點）：
可以看到出現(xiàn)了些不應該出現(xiàn)的頻率，也驗證了我們的推論，其代碼如下：
Fs=512;%采樣頻率
Ts=1/Fs;%采樣周期
N=128;%N—DFT
L=128;%原信號序列長度
t=(0:L-1).*Ts;%時域自變量
x=cos(10*2*pi*t);%原信號
subplot(311);
plot(t,x);
title("原信號")
xlabel("t/s")
grid on
Y=fft(x,N);%fft變換
Y=abs(Y)./N;%實際幅值變換
f=(0:N-1)*Fs./N;%實際頻率變換
subplot(312);
stem(f(1:N/2),Y(1:N./2));
title("N-DFT變換幅頻響應單邊")
xlabel("f/Hz")
grid on
f=f-Fs./2;%移位
subplot(313);
stem(f,fftshift(Y));%移位
title("N-DFT變換幅頻響應雙邊")
xlabel("f/Hz")
grid on
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
3、那么我們再來看下增大N點是否可以減小頻譜泄露的影響，考慮Fs=512，L=128，N=1024，滿足上述三個關系吧，
Fs=512>=2fc；LxTs=0.15s>=T（不是T的整數(shù)倍，故有頻率泄露）；N=1024>=L=128；那么我們看下仿真結果是否對得到（這次頻譜用stem繪制更明顯點）：
很明顯可以看到與N=128時，其往所求正確頻率處更靠近集中，故增大N可以減小頻率響應的影響。代碼如下：
Fs=512;%采樣頻率
Ts=1/Fs;%采樣周期
N=1024;%N—DFT
L=128;%原信號序列長度
t=(0:L-1).*Ts;%時域自變量
x=cos(10*2*pi*t);%原信號
subplot(311);
plot(t,x);
title("原信號")
xlabel("t/s")
grid on
Y=fft(x,N);%fft變換
Y=abs(Y)./N;%實際幅值變換
f=(0:N-1)*Fs./N;%實際頻率變換
subplot(312);
stem(f(1:N/2),Y(1:N./2));
title("N-DFT變換幅頻響應單邊")
xlabel("f/Hz")
grid on
f=f-Fs./2;%移位
subplot(313);
stem(f,fftshift(Y));%移位
title("N-DFT變換幅頻響應雙邊")
xlabel("f/Hz")
grid on
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
4、最后看下它可不可ifft（）變回去：
這是對于1的ifft只在后面加上下述代碼：
figuresubplot(211)plot(t,x)title("原信號")xlabel("t/s")grid onY=fft(x);xx=ifft(Y);subplot(212)plot(t,xx)title("ifft原信號")xlabel("t/s")grid on
1
2
3
4
5
6
7
8
9
10
11
12
13
所以ifft變換時，應把一開始的N因子變回去。
對于非周期信號，其原理差不多，關鍵是時域能否盡可能取到一個周期等，或者對于無限長的信號，只能采取截取出其中最主要的信息了。
綜上，差不多就這樣了，可能多少有些問題，歡迎大家指正互相進步。
