4窗函数及频谱分析

窗函数及频谱分析实验目的：1.掌握各类窗函数的时域和频率特性；2.掌握合理运用窗函数分析信号频谱的方法；3.掌握利用DFT分析连续信号频谱的方法；4.掌握谱分析中参数的选取方法。实验原理：一、窗函数分析在确定信号谱分析中，截短无穷长的序列会造成频率泄漏，合理选取窗函数的类型，可以改善泄露现象。1.常用窗函数矩形窗w=boxcar(N)汉明窗w=hamming(N)汉宁窗w=hanning(N)布莱克曼窗w=blackman(N)凯泽窗w=Kaiser(N,beta)例:N=50;w=boxcar(N);W=fft(w,256);subplot(2,1,1);stem([0:N-1],w);subplot(2,1,2);plot([-128:127],abs(fftshift(W)))MATLAB中提供了fft函数，FFT是DFT的快速算法。X=fft(x,n)：补零或截短的n点傅立叶变换；fftshift(x)将fft计算输出的零频移到输出的中心。例：N=50;w=hamming(N);W=fft(w,256);subplot(2,1,1);stem([0:N-1],w);subplot(2,1,2);plot([-128:127],abs(fftshift(W)))例：已知一连续信号为12()cos(2)cos(2)xtftft其中f1=100Hz，f2=120Hz，若以抽样频率fsam=600Hz对该信号进行抽样，试用DFT近似分析其频谱：利用不同宽度N的矩形窗截短该序列，N分别取15，40，80观察不同长度的窗对谱分析结果的影响；利用汉明窗重做（1）。用矩形窗分析：N=input('请输入N的值：');L=512;f1=100;f2=120;fs=600;ws=2*pi*fs;t=(0:N-1)*(1/fs);x=cos(2*pi*f1*t)+cos(2*pi*f2*t);subplot(211);stem(t,x);W=fft(x,L);f=((-L/2:L/2-1)*(2*pi/L)*fs)/(2*pi);%f=((-L/2:L/2-1)*(1/L)*fs);subplot(212);plot(f,abs(fftshift(W)))用汉明窗重做上述谱分析：N=input('请输入N的值：');L=512;f1=100;f2=120;fs=600;ws=2*pi*fs;t=(0:N-1)*(1/fs);x=cos(2*pi*f1*t)+cos(2*pi*f2*t);wh=hamming(N)';x=x.*wh;subplot(211);stem(t,x);W=fft(x,L);f=((-L/2:L/2-1)*(2*pi/L)*fs)/(2*pi);subplot(212);plot(f,abs(fftshift(W)))例：已知连续信号为12()cos(2)0.15cos(2)xtftft,其中f1=100Hz，f2=150Hz，若以抽样频率fsam=600Hz对该信号进行抽样，利用不同宽度N的矩形窗截短该序列，N分别取15，40，80观察不同长度的窗对谱分析结果的影响；用汉明窗重做上述谱分析。用矩形窗：N=input('请输入N的值：');L=512;f1=100;f2=150;fs=600;ws=2*pi*fs;t=(0:N-1)*(1/fs);x=cos(2*pi*f1*t)+0.15*cos(2*pi*f2*t);subplot(211);stem(t,x);W=fft(x,L);f=((-L/2:L/2-1)*(2*pi/L)*fs)/(2*pi);subplot(212);plot(f,abs(fftshift(W)))用汉明窗：N=input('请输入N的值：');L=512;f1=100;f2=150;fs=600;ws=2*pi*fs;t=(0:N-1)*(1/fs);x=cos(2*pi*f1*t)+0.15*cos(2*pi*f2*t);wh=hamming(N)';x=x.*wh;subplot(211);stem(t,x);W=fft(x,L);f=((-L/2:L/2-1)*(2*pi/L)*fs)/(2*pi);subplot(212);plot(f,abs(fftshift(W)))二、利用DFT计算连续非周期信号的频谱。分析步骤：（1）根据时域抽样定理，确定时域抽样间隔T，得到离散序列[]xk；（2）确定信号截短的长度M及窗函数的类型，得到有限长M点离散序列[][][]Mxkxkwk；（3）确定频域抽样点数N，要求N≥M；（4）利用FFT函数进行N点FFT计算得到N点的[]Xm；（5）由[]Xm可得连续信号的频谱()Xj样点的近似值。例：利用DFT近似分析连续信号x(t)=e-tu(t)的频谱，要求频率分辨率为0.1Hz，绘出频谱图并与理论值比较，。信号的频谱为1()1Xjj，幅度谱为21()1Xj此信号的频谱为无限宽且单调递减，当ω=2π×25rad/s，|X(j50π)|=0.0064,已经衰减到很小，因此选fm=25Hz为近似的最高频率，则抽样频率为fs=2fm=50Hz，抽样间隔为T=1/2fm=0.02s。时域信号无限长，因此必须截短到有限长。由/50/0.1500sNff，则截短点数为N=500，抽样持续时间为Tp=NT=1/1/0.110f，FFT点数可选512，采用矩形窗：fsam=input('请输入fsam的值：');Tp=10;N=512;T=1/fsam;t=0:T:Tp;x=exp(-1*t);X=T*fft(x,N);%消除1/T因子的影响w=(-N/2:N/2-1)*(2*pi/N)*fsam;y=1./(j*w+1);%理论频谱值figure;plot(w,abs(fftshift(X)),w,abs(y),'r-.');title('幅度谱');xlabel('w');legend('计算值','理论值');由幅度谱可知，DFT的计算结果与理论值基本吻合。若将抽样频率降低至10Hz、5Hz等，试观察DFT计算结果与理论值的误差大小。实验内容：1.已知信号表达式为12()cos(2)0.2cos(2)xtftft，其中f1=100Hz，f2=120Hz，试确定合适的抽样频率fs并用DFT近似分析其频谱：（1）利用不同宽度N的矩形窗截短该序列，观察不同长度的窗对谱分析结果的影响，是否都能分辨到两个明显谱峰，为什么？（2）利用汉明窗重做（1），比较与（1）中结果的区别，并解释原因。forN=[15,30,45]L=512;f1=100;f2=120;fs=600;ws=2*pi*fs;t=(0:N-1)*(1/fs);x=cos(2*pi*f1*t)+0.15*cos(2*pi*f2*t);figure;subplot(211);stem(t,x);W=fft(x,L);f=((-L/2:L/2-1)*(2*pi/L)*fs)/(2*pi);subplot(212);plot(f,abs(fftshift(W)))endforN=[15,30,45];L=512;f1=100;f2=120;fs=600;ws=2*pi*fs;t=(0:N-1)*(1/fs);x=cos(2*pi*f1*t)+0.15*cos(2*pi*f2*t);wh=hamming(N)';x=x.*wh;figure;subplot(211);stem(t,x);W=fft(x,L);f=((-L/2:L/2-1)*(2*pi/L)*fs)/(2*pi);subplot(212);plot(f,abs(fftshift(W)))end00.0050.010.0150.020.0250.030.0350.040.0450.05-2-1012-300-200-10001002003000510152000.010.020.030.040.050.060.070.08-2-1012-300-200-1000100200300051015202500.0050.010.0150.020.025-1-0.500.51-300-200-10001002003000123400.0050.010.0150.020.0250.030.0350.040.0450.05-1-0.500.51-300-200-10001002003000246800.010.020.030.040.050.060.070.08-1-0.500.511.5-300-200-10001002003000510152.已知连续信号x(t)=e-2tu(t)，选取合适参数利用DFT近似分析其幅度谱，要求频率分辨率为0.05Hz，绘出谱图并与理论值比较，将理论频谱曲线和实际计算频谱曲线绘制在一个坐标系中。（若将抽样频率降低至10Hz、5Hz等，再重做一次，观察最终理论值与计算值间的误差变化。）forfsam=[50,10,5];Tp=20;N=512;T=1/fsam;t=0:T:Tp;x=exp(-2*t);X=T*fft(x,N);%消除1/T因子的影响w=(-N/2:N/2-1)*(2*pi/N)*fsam;y=1./(j*w+2);%理论频谱值figure;plot(w,abs(fftshift(X)),w,abs(y),'r-.');title('幅度谱');xlabel('w');legend('计算值','理论值');end-200-150-100-5005010015020000.10.20.30.40.50.60.7幅度谱w计算值理论值-40-30-20-1001020304000.10.20.30.40.50.60.7幅度谱w计算值理论值-20-15-10-50510152000.10.20.30.40.50.60.7幅度谱w计算值理论值实验报告要求：1.列出本次实验编写的所有MATLAB程序及各项实验结果的图形（打印），对实验结果进行必要的分析说明；2.总结实验体会及实验中存在的问题。