基于加窗FFT的频谱分析测绘1603班裴浩阳.docx
《基于加窗FFT的频谱分析测绘1603班裴浩阳.docx》由会员分享，可在线阅读，更多相关《基于加窗FFT的频谱分析测绘1603班裴浩阳.docx（37页珍藏版）》请在冰豆网上搜索。
基于加窗FFT的频谱分析测绘1603班裴浩阳
《数字信号处理》
题目：
基于加窗FFT的频谱分析
SpectrumAnalysisBasedonWindowedFFT
姓名裴浩阳
学院交通运输工程学院
班级测绘1603班
学号201628010320
指导老师王小华
长沙理工大学
2018年10月
目录
摘要1
ABSTRACT1
引言3
一、窗函数+FFT对信号进行频谱分析4
二、截断时间长度对频谱分析的影响13
三、采样频率对频谱分析的影响29
参考文献46
附录46
摘要
根据傅里叶原理，任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加。
由此，我们可以运用FFT变换（即快速傅里叶变换），将原来难以处理的时域信号转换成易于分析的频域信号，即对信号进行频谱分析。
并且随着计算机和微电子技术的飞速发展，利用FFT变换对信号进行频谱分析已经广泛应用到各个领域，如匹配滤波、雷达干涉信号处理及语音识别等。
本论文主要利用MATLAB计算工具，分别用矩形窗、汉宁窗、哈明窗以及布莱克曼窗+FFT对给定信号进行频谱分析，并分析各窗函数以及在各窗函数下截断时间长度、采样频率对频谱分析的影响。
关键词：
窗函数、FFT变换、频谱分析、截断时间长度、采样频率
ABSTRACT
AccordingtoFourierprinciple,anycontinuousmeasurementsequenceorsignalcanbeexpressedasaninfinitesuperpositionofsinewavesignalswithdifferentfrequencies.Therefore,wecanuseFFTtransform（fastFouriertransform）,transformingthetime-domainsignal,whichisdifficulttoprocessoriginallyintoafrequency-domainsignaleasytoanalyze,thatis,thesignalspectrumanalysis.Withtherapiddevelopmentofcomputerandmicroelectronicstechnology,spectrumanalysisusingFFTtransformhasbeenwidelyusedinvariousfields,suchasmatchedfiltering,radarinterferencesignalprocessingandspeechrecognition.
ThispapermainlyusesMATLABtoanalyzethespectrumofgivensignalwithrectangularwindow,Hanningwindow,HammingwindowandBlackmanwindow+FFT,andanalyzestheinfluenceofwindowfunction,truncationtimeandsamplingfrequencyonspectrumanalysis.
KeyWords:
Windowfunction,FFTtransform,Spectrumanalysis,Truncationtimelength,Samplingfrequency
引言
数字信号处理中通常是取其有限的时间片段进行分析，而不是对无限长的信号进行测量和运算。
具体做法是从信号中截取一个时间片段，然后对信号进行傅里叶变换、相关分析等数学处理。
信号的分析产生了能量泄露，而用FFT算法计算频谱又产生了栅栏效应。
从原理上讲，这两种误差都是不能消除的。
在FFT分析中为了减少或者消除
频谱能量泄露及栅栏效应，可采用不同的截取函数对信号进行截断，阶段函数为窗函数，简称为窗。
不同的窗函数对频谱分析的影响是不一样的，这主要是因为不同的窗函数产生泄露的大小不一样，频率分辨能力也不一样。
信号的加窗处理，重要的问题是在于根据信号的性质和研究目的来选用窗函数。
另，为方便起见，本文所研究的对象信号统一为：
xt=12*sin（w*t+10*pi/180）+6*sin（3*w*t+20*pi/180）+2.5*sin（5*w*t+40*pi/180）+2*sin（7*w*t+60*pi/180）+sin（9*w*t+80*pi/180）+0.5*sin（11*w*t+90*pi/180）;
W=99*pi。
一、窗函数+FFT对信号进行频谱分析
1.用矩形窗+FFT对信号进行频谱分析
1.1MATLAB程序：
clf;
fs=1500;T=1/fs;%fs为采样速率；
Tp=0.1;N=Tp/T;n=0:
N-1;w=99*pi;t=n*T;
xt=12*sin（w*t+10*pi/180）+6*sin（3*w*t+20*pi/180）+2.5*sin（5*w*t+40*pi/180）+2*sin（7*w*t+60*pi/180）+sin（9*w*t+80*pi/180）+0.5*sin（11*w*t+90*pi/180）;
A=1;xt=A*（xt）/max（abs（xt））;
%将函数值xn的值域限定在[-1,1];
win1=xt.*boxcar（N）';
plot（n,win1）;%矩形窗函数
xlabel（'t'）;ylabel（'x（t）'）;
title（'矩形窗函数'）;axis;
figure;
k=1024;f=fs*（0:
k-1）/k;%f为模拟信号；
xk1=fft（win1,k）;%快速傅里叶变换
xk1=A*xk1/max（abs（xk1））;
%将函数值xk1的值域限定在[-1,1];
plot（f,xk1）;
xlabel（'频率（f）'）;ylabel（'幅值（xk1）'）;
title（'矩形窗函数+FFT频谱分析'）;axis;
1.2窗函数图及频谱图
1.2.1窗函数图
1.2.2频谱图
矩形窗函数对频谱分析的影响：
矩形窗进行频谱分析的优点是主瓣比较集中,缺点是旁瓣较高,并有负旁瓣,导致变换中带进了高频干扰和泄漏，甚至出现负谱现象。
频率识别精度最高，幅值识别精度最低。
2.用汉宁窗+FFT对信号进行频谱分析
2.1MATLAB程序：
clf;
fs=1500;T=1/fs;%fs为采样速率；
Tp=0.1;N=Tp/T;n=0:
N-1;w=99*pi;t=n*T;
xt=12*sin（w*t+10*pi/180）+6*sin（3*w*t+20*pi/180）+2.5*sin（5*w*t+40*pi/180）+2*sin（7*w*t+60*pi/180）+sin（9*w*t+80*pi/180）+0.5*sin（11*w*t+90*pi/180）;
A=1;xt=A*（xt）/max（abs（xt））;
%将函数值xn的值域限定在[-1,1];
win2=xt.*hanning（N）';
plot（n,win2）;%汉宁窗函数
xlabel（'t'）;ylabel（'xt'）;title（'汉宁窗函数'）;axis;
figure;
k=1024;f=fs*（0:
k-1）/k;%f为模拟信号；
xk2=fft（win2,k）;%快速傅里叶变换
xk2=A*xk2/max（abs（xk2））;%将函数值xk1的值域限定在[-1,1];
plot（f,xk2）;xlabel（'频率（f）'）;
ylabel（'幅值（xk1）'）;title（'矩形窗函数+FFT频谱分析'）;axis;
2.2窗函数图及频谱图
2.2.1窗函数图
2.2.2频谱图
汉宁窗函数对频谱分析的影响：
汉宁窗进行频谱分析的优点是主瓣加宽并降低，旁瓣显著减小，可减少泄露，优于矩形窗。
缺点是主瓣加宽相当于分析带宽加宽，频率分辨率下降。
它与矩形窗相比，泄漏、波动都减小了,并且选择性也提高。
3.用哈明窗+FFT对信号进行频谱分析
3.1MATLAB程序：
clf;
fs=1500;T=1/fs;%fs为采样速率；
Tp=0.1;N=Tp/T;n=0:
N-1;w=99*pi;t=n*T;
xt=12*sin（w*t+10*pi/180）+6*sin（3*w*t+20*pi/180）+2.5*sin（5*w*t+40*pi/180）+2*sin（7*w*t+60*pi/180）+sin（9*w*t+80*pi/180）+0.5*sin（11*w*t+90*pi/180）;
A=1;xt=A*（xt）/max（abs（xt））;
%将函数值xn的值域限定在[-1,1];
win3=xt.*hamming（N）';
plot（n,win3）;%哈明窗函数
xlabel（'t'）;ylabel（'xt'）;
title（'哈明窗函数'）;axis;
figure;
k=1024;f=fs*（0:
k-1）/k;%f为模拟信号；
xk3=fft（win3,k）;%快速傅里叶变换
xk3=A*xk3/max（abs（xk3））;
%将函数值xk1的值域限定在[-1,1];
plot（f,xk3）;xlabel（'频率（f）'）;ylabel（'幅值（xk1）'）;title（'哈明窗函数+FFT频谱分析'）;axis;
3.2窗函数图及频谱图
3.2.1窗函数图
3.2.2频谱图
哈明窗函数对频谱分析的影响：
哈明窗与汉宁窗都是余弦窗,又称改进的升余弦窗，不过哈明窗能使旁瓣达到更小。
经分析表明,哈明窗的衰减速度比汉宁窗衰减速度慢。
4.用布莱克曼窗+FFT对信号进行频谱分析
4.1MATLAB程序
clf;
fs=1500;T=1/fs;%fs为采样速率；
Tp=0.1;N=Tp/T;n=0:
N-1;w=99*pi;t=n*T;
xt=12*sin（w*t+10*pi/180）+6*sin（3*w*t+20*pi/180）+2.5*sin（5*w*t+40*pi/180）+2*sin（7*w*t+60*pi/180）+sin（9*w*t+80*pi/180）+0.5*sin（11*w*t+90*pi/180）;
A=1;xt=A*（xt）/max（abs（xt））;
%将函数值xn的值域限定在[-1,1];
win4=xt.*blackman（N）';
plot（n,win4）;%布莱克曼窗函数
xlabel（'t'）;ylabel（'xt'）;
title（'布莱克曼窗函数'）;axis;
figure;
k=1024;f=fs*（0:
k-1）/k;%f为模拟信号；
xk4=fft（win4,k）;%快速傅里叶变换
xk4=A*xk4/max（abs（xk4））;
%将函数值xk1的值域限定在[-1,1];
plot（f,xk4）;xlabel（'频率（f）'）;ylabel（'幅值（xk1）'）;title（'布莱克曼窗函数+FFT频谱分析'）;axis;
4.2窗函数图及频谱图
4.2.1窗函数图
4.2.2频谱图
布莱克曼窗函数对频谱分析的影响：
布莱克曼窗主瓣宽，旁瓣比较低，但等效噪声带宽比汉宁窗要大一点，波动却小一点。
频率识别精度最低，但幅值识别精度最高，有更好的选择性。
综上，泄露与窗函数频谱的两侧旁瓣有关，对于窗函数的选用总的原则是，要从保持最大信息和消除旁瓣的综合效果出发来考虑问题，尽可能使窗函数频谱中的主瓣宽度尽量宽，以获得较陡的过滤带；旁瓣衰减应尽量大，以提高阻带的衰减。
在实际问题分析中，我们要学会根据信号的性质和研究目的来选用合适的窗函数。
二、截断时间长度对频谱分析的影响
1.对于矩形窗，在采样频率一定时，增加截断时间长度，分析截断时间长度对频谱分析的影响。
1.1MATLAB程序
clf;
fs=1500;T=1/fs;%fs为采样速率；
Tp1=0.1;Tp2=0.14;Tp3=0.18;
N1=Tp1/T;N2=Tp2/T;N3=Tp3/T;
n1=0:
N1-1;n2=0:
N2-1;n3=0:
N3-1;
w=99*pi;t1=n1*T;t2=n2*T;t3=n3*T;
xt1=12*sin（w*t1+10*pi/180）+6*sin（3*w*t1+20*pi/180）+2.5*sin（5*w*t1+40*pi/180）+2*sin（7*w*t1+60*pi/180）+sin（9*w*t1+80*pi/180）+0.5*sin（11*w*t1+90*pi/180）;
xt2=12*sin（w*t2+10*pi/180）+6*sin（3*w*t2+20*pi/180）+2.5*sin（5*w*t2+40*pi/180）+2*sin（7*w*t2+60*pi/180）+sin（9*w*t2+80*pi/180）+0.5*sin（11*w*t2+90*pi/180）;
xt3=12*sin（w*t3+10*pi/180）+6*sin（3*w*t3+20*pi/180）+2.5*sin（5*w*t3+40*pi/180）+2*sin（7*w*t3+60*pi/180）+sin（9*w*t3+80*pi/180）+0.5*sin（11*w*t3+90*pi/180）;
A=1;xt1=A*（xt1）/max（abs（xt1））;xt2=A*（xt2）/max（abs（xt2））;xt3=A*（xt3）/max（abs（xt3））;
%将函数值xn的值域限定在[-1,1];
win1=xt1.*boxcar（N1）';win2=xt2.*boxcar（N2）';
win3=xt3.*boxcar（N3）';
k=1024;f=fs*（0:
k-1）/k;%f为模拟信号；
xk1=fft（win1,k）;%快速傅里叶变换
xk2=fft（win2,k）;%快速傅里叶变换
xk3=fft（win3,k）;%快速傅里叶变换
xk1=A*xk1/max（abs（xk1））;xk2=A*xk2/max（abs（xk2））;
xk1=A*xk3/max（abs（xk3））;%将函数值xk1的值域限定在[-1,1];
plot（f,xk1）;xlabel（'频率（f）'）;ylabel（'幅值（xk1）'）;
title（'fs=1500;Tp=0.1矩形窗函数+FFT频谱分析'）;axis;
figure;plot（f,xk2）;xlabel（'频率（f）'）;ylabel（'幅值（xk1）'）;
title（'fs=1500;Tp=0.14矩形窗函数+FFT频谱分析'）;axis;
figure;plot（f,xk3）;xlabel（'频率（f）'）;ylabel（'幅值（xk1）'）;
title（'fs=1500;Tp=0.18矩形窗函数+FFT频谱分析'）;axis;
1.2频谱图
1.2.1fs=1500;Tp=0.1s
1.2.2fs=1500;Tp=0.14s
1.2.3fs=1500;Tp=0.18s
2.对于汉宁窗，在采样频率一定时，增加截断时间长度，分析截断时间长度对频谱分析的影响。
2.1MATLAB程序
clf;
fs=1500;T=1/fs;%fs为采样速率；
Tp1=0.1;Tp2=0.14;Tp3=0.18;
N1=Tp1/T;N2=Tp2/T;N3=Tp3/T;
n1=0:
N1-1;n2=0:
N2-1;n3=0:
N3-1;
w=99*pi;t1=n1*T;t2=n2*T;t3=n3*T;
xt1=12*sin（w*t1+10*pi/180）+6*sin（3*w*t1+20*pi/180）+2.5*sin（5*w*t1+40*pi/180）+2*sin（7*w*t1+60*pi/180）+sin（9*w*t1+80*pi/180）+0.5*sin（11*w*t1+90*pi/180）;
xt2=12*sin（w*t2+10*pi/180）+6*sin（3*w*t2+20*pi/180）+2.5*sin（5*w*t2+40*pi/180）+2*sin（7*w*t2+60*pi/180）+sin（9*w*t2+80*pi/180）+0.5*sin（11*w*t2+90*pi/180）;
xt3=12*sin（w*t3+10*pi/180）+6*sin（3*w*t3+20*pi/180）+2.5*sin（5*w*t3+40*pi/180）+2*sin（7*w*t3+60*pi/180）+sin（9*w*t3+80*pi/180）+0.5*sin（11*w*t3+90*pi/180）;
A=1;xt1=A*（xt1）/max（abs（xt1））;
A=1;xt2=A*（xt2）/max（abs（xt2））;
A=1;xt3=A*（xt3）/max（abs（xt3））;
%将函数值xn的值域限定在[-1,1];
win1=xt1.*hanning（N1）';
win2=xt2.*hanning（N2）';win3=xt3.*hanning（N3）';
k=1024;f=fs*（0:
k-1）/k;%f为模拟信号；
xk1=fft（win1,k）;%快速傅里叶变换
xk2=fft（win2,k）;%快速傅里叶变换
xk3=fft（win3,k）;%快速傅里叶变换
xk1=A*xk1/max（abs（xk1））;
xk2=A*xk2/max（abs（xk2））;
xk1=A*xk3/max（abs（xk3））;
%将函数值xk1的值域限定在[-1,1];
plot（f,xk1）;xlabel（'频率（f）'）;ylabel（'幅值（xk1）'）;
title（'fs=1500;Tp=0.1汉宁窗函数+FFT频谱分析'）;axis;
figure;plot（f,xk2）;xlabel（'频率（f）'）;ylabel（'幅值（xk1）'）;
title（'fs=1500;Tp=0.14汉宁窗函数+FFT频谱分析'）;axis;
figure;plot（f,xk3）;
xlabel（'频率（f）'）;ylabel（'幅值（xk1）'）;
title（'fs=1500;Tp=0.18汉宁窗函数+FFT频谱分析'）;axis;
2.2频谱图
2.2.1fs=1500;Tp=0.1s
2.2.2fs=1500;Tp=0.14s
2.2.3fs=1500;Tp=0.18s
3.对于哈明窗，在采样频率一定时，增加截断时间长度，分析截断时间长度对频谱分析的影响。
3.1MATLAB程序
clf;
fs=1500;T=1/fs;%fs为采样速率；
Tp1=0.1;Tp2=0.14;Tp3=0.18;
N1=Tp1/T;N2=Tp2/T;N3=Tp3/T;
n1=0:
N1-1;n2=0:
N2-1;n3=0:
N3-1;
w=99*pi;t1=n1*T;t2=n2*T;t3=n3*T;
xt1=12*sin（w*t1+10*pi/180）+6*sin（3*w*t1+20*pi/180）+2.5*sin（5*w*t1+40*pi/180）+2*sin（7*w*t1+60*pi/180）+sin（9*w*t1+80*pi/180）+0.5*sin（11*w*t1+90*pi/180）;
xt2=12*sin（w*t2+10*pi/180）+6*sin（3*w*t2+20*pi/180）+2.5*sin（5*w*t2+40*pi/180）+2*sin（7*w*t2+60*pi/180）+sin（9*w*t2+80*pi/180）+0.5*sin（11*w*t2+90*pi/180）;
xt3=12*sin（w*t3+10*pi/180）+6*sin（3*w*t3+20*pi/180）+2.5*sin（5*w*t3+40*pi/180）+2*sin（7*w*t3+60*pi/180）+sin（9*w*t3+80*pi/180）+0.5*sin（11*w*t3+90*pi/180）;
A=1;xt1=A*（xt1）/max（abs（xt1））;
A=1;xt2=A*（xt2）/max（abs（xt2））;
A=1;xt3=A*（xt3）/max（abs（xt3））;
%将函数值xn的值域限定在[-1,1];
win1=xt1.*hamming（N1）';
win2=xt2.*hamming（N2）';win3=xt3.*hamming（N3）';
k=1024;f=fs*（0:
k-1）/k;%f为模拟信号；
xk1=fft（win1,k）;%快速傅里叶变换
xk2=fft（win2,k）;%快速傅里叶变换
xk3=fft（win3,k）;%快速傅里叶变换
xk1=A*xk1/max（abs（xk1））;
xk2=A*xk2/max（abs（xk2））;
xk1=A*xk3/max（abs（xk3））;
%将函数值xk1的值域限定在[-1,1];
plot（f,xk1）;xlabel（'频率（f）'）;ylabel（'幅值（xk1）'）;
title（'fs=1500;Tp=0.1哈明窗函数+FFT频谱分析'）;axis;
figure;plot（f,xk2）;xlabel（'频率（f）'）;ylabel（'幅值（xk1）'）;
title（'fs=1500;Tp=0.14哈明窗函数+FFT频谱分析'）;axis;
figure;plot（f,xk3）;
xlabel（'频率（f）'）;ylabel（'幅值（xk1）'）;
title（'fs=1500;Tp=0.18哈明窗函数+FFT频谱分析'）;axis;
3.2频谱图
3.2.1fs=1500;Tp=0.1s
2.2.2fs=1500;Tp=0.14s
3.2.3fs=1500;Tp=0.18s
4.对于布莱克曼窗，在采样频率一定时，增加截断时间长度，分析截断时间长度对频谱分析的影响。
4.1MATLAB程序
clf;
fs=1500;T=1/fs;%fs为采样速率；
Tp1=0.1;Tp2=0.14;Tp3=0.18;
N1=Tp1/T;N2=Tp2/T;N3=Tp3/T;
n1=0:
N1-1;n2=0:
N2-1;n3=0:
N3-1;
w=99*pi;t1=n1*T;t2=n2*T;t3=n3*T;
xt1=12*sin（w*t1+10*pi/180）+6*sin（3*w*t1+20*pi/180）+2.5*sin（5*w*t1+40*pi/180）+2*sin（7*w*t1+60*pi/180）+sin（9*w*t1+80*pi/180）+0.5*sin（11*w*t1+90*pi/180）;
xt2=12*sin（w*t2+10*pi/1
升级会员 