DSP实验7

实验七：窗函数及其在谱分析中的作用
实验目的：
在理论学习的基础上，掌握不同窗函数的性质、特点，并通过实验认识它们在克服FFT 频谱分析的能量泄漏和栅栏效应误差中的作用，以便在实际工作中能根据具体情况正确选用窗函数。

实验任务：
1.执行下面例程，分析不同窗函数的特点并比较结果。

2.编程实现汉宁窗。

演示其时域和频域波形。

信号的截断产生了能量泄漏，而用FFT算法计算频谱又产生了栅栏效应，从原理上讲这两种误差都是不能消除的，但是我们可以通过选择不同的窗函数对它们的影响进行抑制。

下图是几种常用的窗函数的时域和频域波形，其中矩形窗主瓣窄，旁瓣大，频率识别精度最高，幅值识别精度最低；布莱克曼窗主瓣宽，旁瓣小，频率识别精度最低，但幅值识别精度最高。

下面是四个常用窗函数的示例程序，可在Matlab下执行，注意比较它们的特点。

1.矩形窗：
wp=0.2*pi;ws=0.3*pi;rp=.25;as=50;
delta1=(10^(rp/20)-1)/(10^(rp/20)+1);
delta2=(1+delta1)*(10^(-as/20));
deltah=max(delta1,delta2);
deltal=min(delta1,delta2);
weights=[delta2/delta1 1];
deltaf=(ws-wp)/(2*pi);
M=ceil((-20*log10(sqrt(delta1*delta2))-13)/(14.6*deltaf)+1); M=M+4;f=[0 wp/pi ws/pi 1];m=[1 1 0 0];
h=remez(M-1,f,m,weights);
delta_w=2*pi/1000;wsi=ws/delta_w+1;wpi=wp/delta_w;
asd=-max(db(wsi:501));
figure(4);
subplot(211);stem([0:M-1],h);title('actual impulse response') axis([0 M-1 -0.1 0.3]),xlabel('n');ylabel('h(n)')
set(gca,'xtick',[0,M-1])
set(gca,'ytick',[-0.1 0 0.1 0.2 0.3]),grid
subplot(212);plot(w/pi,db),title('magnitude response in db'), axis([0 1 -80 10]),xlabel('frequency in pi units')
ylabel('decibels')
set(gca,'xtick',[0 .2 .3 1])
set(gca,'ytick',[-50,0])
%set(gca,'YTickLabels',['50';'0'])
grid
2.汉明窗
wp=0.2*pi;ws=0.3*pi;tr_width=ws-wp;
m=ceil(8*pi/tr_width)+1;n=[0:m-1];
wc=(ws+wp)/2;hd=ideal_lp(wc,m);
w_ham=(hamming(m))';
h=hd.*w_ham;
delta_w=2*pi/1000;
rp=-(min(db(1:1:wp/delta_w+1)));
as=-round(max(db(ws/delta_w+1:1:501)));
figure(1)
subplot(221);stem(n,hd);title('ideal impuse response')
axis([0 m-1 -.1 .3]);ylabel('hd(n)')
subplot(222);stem(n,w_ham);title('hamming window')
axis([0 m-1 0 1.1]);ylabel('w(n)')
subplot(223);stem(n,h);title('actual impulse response')
axis([0 m-1 -0.1 .3]);ylabel('h(n)')
subplot(224);plot(w/pi,db);title('magnitude response in db');grid axis([0 1 -100 10]);ylabel('decibels')
set(gca,'xtickmode','manual','xtick',[0,.2,.3,1])
set(gca,'ytickmode','manual','ytick',[-50,0])
%set(gca,'yticklabelmode','manual','yticklabels',['50';'0'])
3.布莱克曼窗
ws1=.2*pi;wp1=.35*pi;wp2=.65*pi;ws2=.8*pi;
as=60;
tr_width=min((wp1-ws1),(ws2-wp2));
m=ceil(12*pi/tr_width)+1;
n=[0:m-1];
wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;
hd=ideal_lp(wc2,m)-ideal_lp(wc1,m);
w_bla=(blackman(m))'; h=hd.*w_bla;
delta_w=2*pi/1000;
rp=-min(db(wp1/delta_w+1:wp2/delta_w));
as=-round(max(db(ws2/delta_w+1:501)));
figure(2)
subplot(221);stem(n,hd);title('ideal impulse response')
axis([0 m-1 -0.4 0.5]);ylabel('hd(n)')
subplot(222);stem(n,w_bla);title('blackman window')
axis([0 m-1 0 1.1]);ylabel('w(n)')
subplot(223);stem(n,h);title('actual impulse response')
axis([0 m-1 -0.4 0.5]);ylabel('h(n)')
subplot(224);plot(w/pi,db);title('magnitude response in db');grid ylabel('decibels'),axis([0 1 -150 10])
set(gca,'xtickmode','manual','xtick',[0,.2,.35,.65,.8,1])
set(gca,'ytickmode','manual','ytick',[-60,0])
%set(gca,'yticklabelmode','manual','yticklabels',['60';'0'])
4.凯泽窗
m=45;as=60;n=[0:m-1];
beta=0.1102*(as-8.7);
w_kai=(kaiser(m,beta))';
wc1=pi/3;wc2=2*pi/3;
hd=ideal_lp(wc1,m)+ideal_lp(pi,m)-ideal_lp(wc2,m);
h=hd.*w_kai;
figure(3)
subplot(2，2，1);stem(n,hd);title('ideal impulse response') axis([-1 m -0.2 0.8]);xlabel('n');ylabel('hd(n)')
subplot(2，2，2);stem(n,w_kai);title('kaiser window')
axis([-1 m 0 1.1]);xlabel('n');ylabel('w(n)')
subplot(2，2，3);stem(n,h);title('actual impulse response') axis([-1 m -0.2 0.8]);xlabel('n');ylabel('h(n)')
subplot(2，2，4);plot(w/pi,db)
title('magnitude response in db');grid;
xlabel('frequency in pi units');ylabel('decibels') axis([0 1 -80 10]),
set(gca,'xtickmode','manual','xtick',[0,1/3,2/3,1]) set(gca,'ytickmode','manual','ytick',[-60,0])。