FFT詳解
- 頻域模值*2/N=時域幅度
- fft/ifft前后的數(shù)據(jù)會在能量上差
sqrt(n); 可以使用rms()來查看；
pow_pre=rms(data);
pow_post =rms(fft(data));
ratio=pow_post/pow_pre;
- fft的頻率分辨率為
fs/N - 第n個點的頻率為
Fn=(n-1)fs/N - fft運(yùn)算本質(zhì)上是一個矩陣乘法,下面兩個運(yùn)算
a=b
a=fft(data)
b=data*dftmtx(length(data))
clear all;
N=2048; %fft點數(shù)
fs=512; %采樣率
ts=1/fs;
t=0:ts:(N-1)*ts;
d=10*cos(2*pi*20*t)+2*cos(2*pi*150*t);
df=fftshift(fft(d)).*2/N; %fft前后的幅度變化
f_x = [-N/2:N/2-1]*fs/N;
plot(f_x,df);
xlabel('frequency (Hz)')
ylabel('Magnitude')
1.fft點數(shù)不同對FFT的影響
FFT的結(jié)果在-fs/2：fs/2-1之間奸攻，共有N個點褒墨。FFT的點數(shù)代表頻率分辨率；頻率分辨率是f0=fs/N。
不同的fs之間fft的區(qū)別；fs的增大意味著分析帶寬的增大；同時會降低頻率分辨率轿塔；
頻率分辨率f0=fs/N；可見降低fs或增大N都可以使得f0減小，即提高頻率分辨率催训；
時域補(bǔ)零相當(dāng)于頻域插值
close all; clc; clear all;
%% 128point fft vs 128point with zeros padding
fs=1e3;
N=128;
t=[0:N-1]/fs;
f1=200;f2=100;
a1=0.5;a2=1;
x=a1*cos(2*pi*f1*t)+a2*sin(2*pi*f2*t);
x_len=length(x);
W=dftmtx(N);
x_fft = x*W;
x_max = max(x_fft);
f=(-1/2*N:1/2*N-1)*fs/N;
figure
subplot(2,2,1)
plot(f,fftshift(abs(x_fft)));
title('128point')
subplot(222)
t_less=[0:N-29]/fs;
x_less=a1*cos(2*pi*f1*t_less)+a2*sin(2*pi*f2*t_less);
x_p = [x,zeros(1,128)];
x_p_fft = fft(x_p);
f=(-1/2*256:1/2*256-1)*fs/256;
plot(f,abs(fftshift(x_p_fft)))'
title('zeros padding fft256');
%% 128fft vs 512fft
fs=1e3;
N=128;
t=[0:N-1]/fs;
f1=200;f2=100;
a1=0.5;a2=1;
x=a1*cos(2*pi*f1*t)+a2*sin(2*pi*f2*t);
x_len=length(x);
W=dftmtx(N);
x_fft = x*W;
x_max = max(x_fft);
f=(-1/2*N:1/2*N-1)*fs/N;
subplot(223)
plot(f,fftshift(abs(x_fft)));
title('128point')
N=512;
t=[0:N-1]/fs;
x=a1*cos(2*pi*f1*t)+a2*sin(2*pi*f2*t);
W=dftmtx(N);
x_fft = x*W;
f=(-1/2*N:1/2*N-1)*fs/N;
subplot(224)
plot(f,fftshift(abs(x_fft)));
title('512point');
%% 128point fft vs 128point with zeros padding
fs=1e3;
N=512;
t=[0:N-1]/fs;
f1=200;f2=100;
a1=0.5;a2=1;
x=a1*cos(2*pi*f1*t)+a2*sin(2*pi*f2*t);
x_len=length(x);
W=dftmtx(N);
x_fft = x*W;
x_max = max(x_fft);
f=(-1/2*N:1/2*N-1)*fs/N;
figure
subplot(1,2,1)
plot(f,fftshift(abs(x_fft)));
title('512point&fs=1e3')
subplot(122)
fs=1e4;
t =[0:N-1]/fs;
x=a1*cos(2*pi*f1*t)+a2*sin(2*pi*f2*t);
x_fft = x*W;
f=(-1/2*N:1/2*N-1)*fs/N;
plot(f,abs(fftshift(x_fft)))'
title('512point&fs=1e4')
2. 采樣率的影響：
因為fft的點數(shù)要是2n;如果fs是信號頻率的2n倍時洽议；所取做fft的點數(shù)n采樣不是整周期采樣；因為fft是把階段信號做周期化處理進(jìn)行的漫拭；如果不是整周期采樣亚兄，則周期化時，首尾采樣點跳變采驻，會引入高頻信號审胚；
fs=1e3;
N=128;
t=[0:N-1]/fs;
f1=200;f2=100;
a1=0.5;a2=1;
x=a1*cos(2*pi*f1*t)+a2*sin(2*pi*f2*t);
x_len=length(x);
W=dftmtx(N);
x_fft = x*W;
x_max = max(x_fft);
f=(-1/2*N:1/2*N-1)*fs/N;
figure
subplot(2,2,1)
plot(f,fftshift(abs(x_fft)));
title('16point & fs=1e3')
subplot(222)
fs=800;
N=128;
t=[0:N-1]/fs;
x=a1*cos(2*pi*f1*t)+a2*sin(2*pi*f2*t);
x_fft=fft(x);
f=(-1/2*N:1/2*N-1)*fs/N;
plot(f,abs(fftshift(x_fft)));
title('16point & fs=800');
subplot(223)
fs=1200;
N=128;
t=[0:N-1]/fs;
x=a1*cos(2*pi*f1*t)+a2*sin(2*pi*f2*t);
x_fft=fft(x);
f=(-1/2*N:1/2*N-1)*fs/N;
plot(f,abs(fftshift(x_fft)));
title('16point & fs=1200');
subplot(224)
fs=2000;
N=128;
t=[0:N-1]/fs;
x=a1*cos(2*pi*f1*t)+a2*sin(2*pi*f2*t);
x_fft=fft(x);
f=(-1/2*N:1/2*N-1)*fs/N;
plot(f,abs(fftshift(x_fft)));
title('16point & fs=2000');
3. 窗函數(shù)的影響：
實際情況下，我們得到的信號都是有限長的礼旅，即對原始序列作加窗處理使其成為有限長膳叨，時域的乘積對應(yīng)頻域卷積，造成頻率的泄露痘系。減小泄露的防范可以取更慘的數(shù)據(jù)菲嘴，缺點是運(yùn)算量更大；也可以選擇窗的形狀汰翠，使得窗譜的旁瓣能量更小龄坪。（濾波器是時域卷積相當(dāng)于頻域相乘，加窗和過濾波器運(yùn)算不相同）复唤。
在對信號做FFT分析時健田，如果采樣頻率固定不變，由于被采樣信號自身頻率的微小變化以及干擾因素的影響佛纫，就會使數(shù)據(jù)窗記錄的不是整數(shù)個周期妓局。從時域來說，這種情況在信號的周期延拓時就會導(dǎo)致其邊界點不連續(xù)呈宇，使信號附加了高頻分量; 從頻域來說好爬，由于FFT算法只是對有限長度的信號進(jìn)行變換，有限長度信號在時域相當(dāng)于無限長信號和矩形窗的乘積甥啄，也就是將這個無限長信號截短存炮，對應(yīng)頻域的傅里葉變換是實際信號傅里葉變換與矩形窗傅里葉變換的卷積。
增加采樣長度可以分析出更多頻率的信號型豁，可以減少頻譜泄露，不過增加采樣長度必然會對數(shù)據(jù)處理的實時性造成影響尚蝌！理想的窗函數(shù)是主瓣很窄迎变，旁瓣衰減很快，矩形窗的主瓣很窄飘言，但是旁瓣衰減卻很慢衣形，hanning窗、hamming窗、blackman窗等的旁瓣衰減有了明顯的改進(jìn)谆吴，但是主瓣卻寬了很多倒源，大概是矩形窗主瓣的二倍，blackman窗的主瓣還要寬句狼，這就造成了信號頻譜的頻率識別率很低笋熬！
加窗的區(qū)別：1是沒加窗相當(dāng)于矩形窗，2是hamming腻菇；3是加矩形窗胳螟；
%% 128point fft windows
fs=1e3;
N=128;
t=[0:N-1]/fs;
f1=200;f2=100;
a1=0.5;a2=1;
x=a1*cos(2*pi*f1*t)+a2*sin(2*pi*f2*t);
x_len=length(x);
W=dftmtx(N);
x_fft = x*W;
x_max = max(x_fft);
f=(-1/2*N:1/2*N-1)*fs/N;
figure
subplot(2,2,1)
plot(f,fftshift(abs(x_fft)));
title('512point&fs=1e3')
subplot(222)
win=hamming(N);
x_fft = x.*win'*W;
plot(f,fftshift(abs(x_fft)));
title('hamming')
subplot(223)
rwin = rectwin(N);
x_fft = x.*rwin'*W;
plot(f,fftshift(abs(x_fft)));
title('rectwin')
參考資料：
http://blog.jobbole.com/70549/
http://www.cnblogs.com/v-July-v/archive/2011/02/20/1983676.html
