窗函数性能分析——MATLAB

一、实验目的

1、掌握Matlah由各种窗数序列的生成方法;
2、掌握窗函数序列频率特性的计算与画图方法;
3、掌握窗函数的相关参数对窗函数频域性能的影响;
4、了解混合窗函数的定义、生成方法和频域性能特点;
5、对各种常用的窗函数序列的时域及频域性能特点有整体的认识;
6、了解窗函数序列实现频率响应的线性相位特性需要满足的基本条件。

二、实验原理

1、常用窗函数序列

加窗处理是信号处理中常用的一种处理方式，在数字滤波器的设计、信号的离散谱分析中都有很广泛的应用。常用的窗函数序列主要有:

(1）矩形窗函数序列

矩形窗序列是最基本形式的一种窗函数序列，其时域表示如
公式所示:

(2）汉宁窗函数序列

汉宁窗序列也称为升余弦窗，在矩形窗序列的基础上，增加一个余弦项，其时域表示如公式所示

(3）海明窗函数序列

海明窗序列和汉宁窗序列类似，只是窗函数的两部分的系数有所不同,其时域表示如公式所示:

(4）布莱克曼窗序列

布莱克曼窗序列进一步增加一个余弦项,其时域表示如公式所示:

(5）凯瑟窗序列

凯瑟窗序列和前面4种类型的窗函数序列在原理上有所不同,
其建立在第一类修正的零阶贝塞尔函数基础上,其时域表示如公式所示:

公式中，I0(g)是第一类修正的零阶贝塞尔函数，β为调整参数,用以调整窗函数序列的性能
在 Matlab中,各种窗函数序列的生成均有相对应的函数，对应上述的各窗函数序列的分别为;
(1）矩形窗序列:boxcar()
(2)汉宁窗序列:hanning()
(3）海明窗序列:hamming()(汉明窗)
(4）布莱克曼窗序列:blackman()
(5）凯瑟窗序列:kaiser()
这几个窗函数的调用格式比较简单,其中,前4个函数只要给出窗函数的长度N即可,如:

Window=hamming(N);

该函数产生一个长度为N的海明窗序列,结果保存在输出向量Window中。
kaiser()函数在使用中不仅要给出窗序列的长度N,还需给出窗函数的参数β，即其调用格式为:

Window=kaiser(N,belta);

其中,输入参数belta即为窗函数的参数β，一般在4~8之间取值,用以调整凯瑟窗序列的性能。

2、窗函数的频域特性

将各种窗函数的时域表示式进行离散时间傅里叶变换,即可得到窗函数的频域表示。各种窗函数的幅频特性一般都具有下图的形式:

衡量窗函数的性能,通常通过两个指标:

（1）主瓣宽度。

即窗函数幅频特性曲线第一至第二过零点之间的频带宽度。在滤波器设计中,主瓣越宽，则滤波器过渡带越宽;在频谱分析中，主瓣越宽，则频率分辨力越差。所以在实际应用中,总是希望窗函数有尽量窄的主瓣宽度。

(2）副瓣衰减。

一般用第一副瓣相对于主瓣的衰减进行表示。在滤波器设计中,副瓣衰减越大，则滤波器幅频响应在通带之内的平稳程度和阻带衰减特性就越好;在频谱分析中，副瓣衰减越大,则频谱泄漏效应就会减弱。所以在实际应用中，总是希望窗函数有尽量大的副瓣衰减。窗函数的主全宽度和副瓣衰减性能之间往往是相互矛盾的,提高了一个性能,另外一个性能就会有所下降。
在Matlab中,窗函数序列频率特性的计算可由 freqz）函数来实现。以下以汉宁窗序列为例，给出窗函数序列的频率特性计算及画图的仿真程序实现方法,如下:

N=10;%设置窗函数的宽度成汉宁窗
wn=hanning(N);%生成频率特性
[H,m] = freqz(wn, [1],1024,'whole');%计算窗函数的算幅度特性
mag = abs(H);
db =20*log10((mag"+eps)/max(mag));%将幅度特性转换为对数形式算相位特性
pha=angle(H);
subplot(121)%以下代码为绘制幅度特性曲线plot(m/pi,db); axis([01 -100 0]);
xlabel('w/pi');ylabel('dB');title('幅度特性); grid on;
subplot(122)%以下代码为绘制相位特性曲线plot(m/pi,pha); axis([0 1 -pi pi]);
xlabel ('w/pi');ylabel('rad');title('相位特性'); grid on

各种窗函数序列的相位特性一般都是线性的（或分段线性）,即为频率的线性函数,也称为线性相位特性。但是,要保证窗函数的线性相位特性,需要保证窗函数序列是中心对称的,即需要满足:

3、混合窗函数序列简介

混合窗函数序列是目前在信号处理、滤波器设计等领域经常提及的一种窗函数序列,其幅度特性较常用的窗函数序列更为理想一些。滤合窗承数序列的定义如下所示:

公式中，窗函数的前半部分为矩形窗利叠加，而后半部分则是标准的海明窗的后一半。α、β和γ是可调整系数,会影响窗函数的主那宽度和副瓣衰减, 这三个参数在实际应用中可根据实际需求进行调整，但需满足α+β+γ=1。由混称的,因而其相合窗函数的定义也可看出混合窗函数是非中心对位特性是非线性的,从本质上讲是一种以牺牲相位特性来提升幅度特性的窗函数序列

三、实验步骤、数据记录及处理

本实验以Matlab信号处理工具箱所提供的窗函数为基础，结合频率特性分析,对常用窗函数的时域和频域特性进行分析和对验内容和实验步比,并对相关的结论与参数进行验证。具体的步骤如下:
1、设置窗函数的长度N=15,分别生成矩形窗、汉宁窗、海明窗和布莱克曼窗序列，绘制各窗函数序列的时域图形,观察各种窗函数的时域性能。
2、利用频率响应计算的函数或傅立叶变换函数,分别计算实验内容1中所生成的各窗函数序列的频率特性,绘制幅频和相频时对各种窗函数响应曲线，验证理论课程中的相关性能指标，同的频域性能进行对比。
3、设置窗函数的长度N=33,重复上述过程, 分析窗序列的长度对窗函数频域性能的影响。
4、生成N=33，β分别为B=4,5,6,7,8的凯瑟窗,分别计算并绘制幅频特性曲线，分析调整参数β对窗函数频域性能的影响。
5、设置窗序列长度N=31,生成汉宁窗。设置N=31,a=0.42,，生成混合窗序列。画出这两个窗序列的时域图β=0.4，和y=0.18形，并进行对比。
6、计算并绘制5中生成的两个窗函数的幅度特性和相位特性,对两种窗函数的幅度性能进行对比。同时,观察两个窗函数的相位特性,验证窗函数取得线性相位特性需满足的条件。
实验程序：

clc;clear all;
N=15;%设置窗函数序列长度；
%矩形窗序列
window1=boxcar(N);%生成窗函数序列；
[H1,m1]=freqz(window1,[1],1024,'whole');%计算函数的频率特性
mag1=abs(H1);%计算函数的幅度特性；
pha1=angle(H1);%计算函数的相位特性；
db1=20*log10((mag1+eps)/max(mag1));%将幅度特性转换为对数形式
%绘制矩形窗频域波形
figure('name','矩形窗频域波形');
subplot(311);stem(window1);grid on;xlabel('n');ylabel('x(n)');title('N=15矩形窗函数时域图形');
subplot(312);plot(m1/(2*pi),db1);grid on;xlabel('w/(2*pi)');ylabel('ab');title('N=15矩形窗幅频响应');
subplot(313);plot(m1/(2*pi),pha1);grid on;xlabel('w/(2*pi)');ylabel('rad');title('N=15矩形窗相频响应');

%汉宁窗序列
window2=hanning(N);%生成窗函数序列；
[H2,m2]=freqz(window2,[1],1024,'whole');%计算函数的频率特性
mag2=abs(H2);%计算函数的幅度特性；
pha2=angle(H2);%计算函数的相位特性；
db2=20*log10((mag2+eps)/max(mag2));%将幅度特性转换为对数形式
%绘制汉宁窗频域波形
figure('name','汉宁窗频域波形');
subplot(311);stem(window2);grid on;xlabel('n');ylabel('x(n)');title('N=15汉宁窗函数时域图形');
subplot(312);plot(m2/(2*pi),db1);grid on;xlabel('w/(2*pi)');ylabel('ab');title('N=15汉宁窗幅频响应');
subplot(313);plot(m2/(2*pi),pha1);grid on;xlabel('w/(2*pi)');ylabel('rad');title('N=15汉宁窗相频响应');

%汉明窗序列
window3=hamming(N);%生成窗函数序列；
[H3,m3]=freqz(window3,[1],1024,'whole');%计算函数的频率特性
mag3=abs(H3);%计算函数的幅度特性；
pha3=angle(H3);%计算函数的相位特性；
db3=20*log10((mag3+eps)/max(mag3));%将幅度特性转换为对数形式
%绘制汉明窗频域波形
figure('name','汉明窗频域波形');
subplot(311);stem(window3);grid on;xlabel('n');ylabel('x(n)');title('N=15汉明窗函数时域图形');
subplot(312);plot(m3/(2*pi),db1);grid on;xlabel('w/(2*pi)');ylabel('ab');title('N=15汉明窗幅频响应');
subplot(313);plot(m3/(2*pi),pha1);grid on;xlabel('w/(2*pi)');ylabel('rad');title('N=15汉明窗相频响应');

%布莱克曼窗序列
window4=blackman(N);%生成窗函数序列；
[H4,m4]=freqz(window4,[1],1024,'whole');%计算函数的频率特性
mag4=abs(H4);%计算函数的幅度特性；
pha4=angle(H4);%计算函数的相位特性；
db4=20*log10((mag4+eps)/max(mag4));%将幅度特性转换为对数形式
%绘制布莱克曼频域波形
figure('name','布莱克曼频域波形');
subplot(311);stem(window4);grid on;xlabel('n');ylabel('x(n)');title('N=15布莱克曼窗函数时域图形');
subplot(312);plot(m4/(2*pi),db1);grid on;xlabel('w/(2*pi)');ylabel('ab');title('N=15布莱克曼窗幅频响应');
subplot(313);plot(m4/(2*pi),pha1);grid on;xlabel('w/(2*pi)');ylabel('rad');title('N=15布莱克曼窗相频响应');

clc;clear all;
N=33;%设置窗函数序列长度；
%矩形窗函数
window1=boxcar(N);%生成窗函数序列；
[H1,m1]=freqz(window1,[1],1024,'whole');%计算函数的频率特性
mag1=abs(H1);%计算函数的幅度特性；
pha1=angle(H1);%计算函数的相位特性；
db1=20*log10((mag1+eps)/max(mag1));%将幅度特性转换为对数形式
%绘制矩形窗频域波形
figure('name','矩形窗频域波形');
subplot(311);stem(window1);grid on;xlabel('n');ylabel('x(n)');title('N=33矩形窗时域图形');
subplot(312);plot(m1/(2*pi),pha1);grid on;xlabel('w/2pi');ylabel('ψ(w)');title('N=33矩形窗的相频特性');
subplot(313);plot(m1/(2*pi),db1);grid on;xlabel('w/2pi');ylabel('db');title('N=33矩形窗的幅度特性');

%汉宁窗函数
window2=hanning(N);%生成窗函数序列；
[H2,m2]=freqz(window2,[1],1024,'whole');%计算函数的频率特性
mag2=abs(H2);%计算函数的幅度特性；
pha2=angle(H2);%计算函数的相位特性；
db2=20*log10((mag2+eps)/max(mag2));%将幅度特性转换为对数形式
%绘制汉宁窗频域波形
figure('name','汉宁窗频域波形');
subplot(311);stem(window2);grid on;xlabel('n');ylabel('x(n)');title('N=33的汉宁窗函数时域图形');
subplot(312);plot(m2/(2*pi),pha2);grid on;xlabel('w/2pi');ylabel('ψ(w)');title('N=33的汉宁窗函数的相频特性');
subplot(313);plot(m2/(2*pi),db2);grid on;xlabel('w/2pi');ylabel('db');title('N=33的汉宁窗函数的幅度特性');

%汉明窗函数
window3=hamming(N);%生成窗函数序列；
[H3,m3]=freqz(window3,[1],1024,'whole');%计算函数的频率特性
mag3=abs(H3);%计算函数的幅度特性；
pha3=angle(H3);%计算函数的相位特性；
db3=20*log10((mag3+eps)/max(mag3));%将幅度特性转换为对数形式
%绘制汉明窗频域波形
figure('name','汉明窗频域波形');
subplot(311);stem(window3);grid on;xlabel('n');ylabel('x(n)');title('N=33的汉明窗函数时域图形');
subplot(312);plot(m3/(2*pi),pha3);grid on;xlabel('w/2pi');ylabel('ψ(w)');title('N=33的汉明窗函数的相频特性');
subplot(313);plot(m3/(2*pi),db3);grid on;xlabel('w/2pi');ylabel('db');title('N=33的汉明窗函数的幅度特性');

%布莱克曼窗函数
window4=blackman(N);%生成窗函数序列；
[H4,m4]=freqz(window4,[1],1024,'whole');%计算函数的频率特性
mag4=abs(H4);%计算函数的幅度特性；
pha4=angle(H4);%计算函数的相位特性；
db4=20*log10((mag4+eps)/max(mag4));%将幅度特性转换为对数形式
%绘制布莱克曼频域波形
figure('name','布莱克曼频域波形');
subplot(311);stem(window4);grid on;xlabel('n');ylabel('x(n)');title('N=33的布莱克曼窗函数时域图形');
subplot(312);plot(m4/(2*pi),pha4);grid on;xlabel('w/2pi');ylabel('ψ(w)');title('N=33的布莱克曼窗函数的相频特性');
subplot(313);plot(m4/(2*pi),db4);grid on;xlabel('w/2pi');ylabel('db');title('N=33的布莱克曼窗函数的幅度特性');

clc;clear all;
N=33;%设置窗函数序列长度；
belta1=4;
window1=kaiser(N,belta1);%凯瑟窗函数
[H1,m1]=freqz(window1,[1],1024,'whole');%计算函数的频率特性
mag1=abs(H1);%计算函数的幅度特性；
pha1=angle(H1);%计算函数的相位特性；
db1=20*log10((mag1+eps)/max(mag1));%将幅度特性转换为对数形式
figure('name','=4时凯瑟窗函数频域波形');
subplot(311);stem(window1);grid on;xlabel('n');ylabel('x(n)');title('β=4的凯瑟窗函数时域图形');
subplot(312);plot(m1/(2*pi),pha1);grid on;xlabel('w/2pi');ylabel('ψ(w)');title('β=4的凯瑟窗函数的相频特性');
subplot(313);plot(m1/(2*pi),db1);grid on;xlabel('w/2pi');ylabel('db');title('β=4的凯瑟窗函数的幅度特性');

belta2=5;
window2=kaiser(N,belta2);%凯瑟窗函数
[H2,m2]=freqz(window2,[1],1024,'whole');%计算函数的频率特性
mag2=abs(H2);%计算函数的幅度特性；
pha2=angle(H2);%计算函数的相位特性；
db2=20*log10((mag2+eps)/max(mag2));%将幅度特性转换为对数形式
figure('name','=5时凯瑟窗函数频域波形');
subplot(311);stem(window2);grid on;xlabel('n');ylabel('x(n)');title('β=5的凯瑟窗函数时域图形');
subplot(312);plot(m2/(2*pi),pha2);grid on;xlabel('w/2pi');ylabel('ψ(w)');title('β=5的凯瑟窗函数的相频特性');
subplot(313);plot(m2/(2*pi),db2);grid on;xlabel('w/2pi');ylabel('db');title('β=5的凯瑟窗函数的幅度特性');

belta3=6;
window3=kaiser(N,belta3);%凯瑟窗函数
[H3,m3]=freqz(window3,[1],1024,'whole');%计算函数的频率特性
mag3=abs(H3);%计算函数的幅度特性；
pha3=angle(H3);%计算函数的相位特性；
db3=20*log10((mag3+eps)/max(mag3));%将幅度特性转换为对数形式
figure('name','=6时凯瑟窗函数频域波形');
subplot(311);stem(window3);grid on;xlabel('n');ylabel('x(n)');title('β=6的凯瑟窗函数时域图形');
subplot(312);plot(m3/(2*pi),pha3);grid on;xlabel('w/2pi');ylabel('ψ(w)');title('β=6的凯瑟窗函数的相频特性');
subplot(313);plot(m3/(2*pi),db3);grid on;xlabel('w/2pi');ylabel('db');title('β=6的凯瑟窗函数的幅度特性');

belta4=7;
window4=kaiser(N,belta4);%凯瑟窗函数
[H4,m4]=freqz(window4,[1],1024,'whole');%计算函数的频率特性
mag4=abs(H4);%计算函数的幅度特性；
pha4=angle(H4);%计算函数的相位特性；
db4=20*log10((mag4+eps)/max(mag4));%将幅度特性转换为对数形式
figure('name','=7时凯瑟窗函数频域波形');
subplot(311);stem(window4);grid on;xlabel('n');ylabel('x(n)');title('β=7的凯瑟窗函数时域图形');
subplot(312);plot(m4/(2*pi),pha4);grid on;xlabel('w/2pi');ylabel('ψ(w)');title('β=7的凯瑟窗函数的相频特性');
subplot(313);plot(m4/(2*pi),db4);grid on;xlabel('w/2pi');ylabel('db');title('β=7的凯瑟窗函数的幅度特性');

belta5=8;
window5=kaiser(N,belta5);%凯瑟窗函数
[H5,m5]=freqz(window5,[1],1024,'whole');%计算函数的频率特性
mag5=abs(H5);%计算函数的幅度特性；
pha5=angle(H5);%计算函数的相位特性；
db5=20*log10((mag5+eps)/max(mag5));%将幅度特性转换为对数形式
figure('name','=8时凯瑟窗函数频域波形');
subplot(311);stem(window5);grid on;xlabel('n');ylabel('x(n)');title('β=8的凯瑟窗函数时域图形');
subplot(312);plot(m5/(2*pi),pha5);grid on;xlabel('w/2pi');ylabel('ψ(w)');title('β=8的凯瑟窗函数的相频特性');
subplot(313);plot(m5/(2*pi),db5);grid on;xlabel('w/2pi');ylabel('db');title('β=8的凯瑟窗函数的幅度特性');

clear all;clc;
N=31;n=0:30;
a=0.42;b=0.4;c=0.18;
x=hanning(N);%汉宁窗函数
for n=0:1:(N-1)/2 %混合窗函数
    w(n+1)=a-b*cos(2*pi*n/(N-1))+c*sin(pi*n/(N-1));
end
for n=(N+1)/2:1:N-1
    w(n+1)=0.54-0.46*cos(2*pi*n/(N-1));
end
%绘制汉宁窗频域波形
figure('name','汉宁窗时域波形');
subplot(211);stem(x);grid on;xlabel('n');ylabel('x(n)');title('N=31的汉宁窗函数时域图形');
subplot(212);stem(w);grid on;xlabel('n');ylabel('w(n)');title('N=31的混合窗函数时域图形');

clear all;clc;
N=31;n=0:30;
a=0.42;b=0.4;c=0.18;
x=hanning(N);%汉宁窗函数
for n=0:1:(N-1)/2 %混合窗函数
    w(n+1)=a-b*cos(2*pi*n/(N-1))+c*sin(pi*n/(N-1));
end
for n=(N+1)/2:1:N-1
    w(n+1)=0.54-0.46*cos(2*pi*n/(N-1));
end
[H1,m1]=freqz(x,[1],1024,'whole');%计算函数的频率特性
mag1=abs(H1);%计算函数的幅度特性；
pha1=angle(H1);%计算函数的相位特性；
db1=20*log10((mag1+eps)/max(mag1));%将幅度特性转换为对数形式
%绘制汉宁窗频域波形
figure('name','汉宁窗频域波形');
subplot(221);plot(m1/(2*pi),pha1);grid on;xlabel('w/2pi');ylabel('ψ(w)');title('N=31的汉宁窗函数的相频特性');
subplot(223);plot(m1/(2*pi),db1);grid on;xlabel('w/2pi');ylabel('db');title('N=31的汉宁窗函数的幅度特性');
[H2,m2]=freqz(w,[1],1024,'whole');%计算函数的频率特性
mag2=abs(H2);%计算函数的幅度特性；
pha2=angle(H2);%计算函数的相位特性；
db2=20*log10((mag2+eps)/max(mag2));%将幅度特性转换为对数形式
subplot(222);plot(m2/(2*pi),pha2);grid on;xlabel('w/2pi');ylabel('ψ(w)');title('N=31的混合窗函数的相频特性');
subplot(224);plot(m2/(2*pi),db2);grid on;xlabel('w/2pi');ylabel('db');title('N=31的混合窗函数的幅度特性');

四、分析

分析略

五、思考题

(1)汉宁窗、海明窗和布莱克曼窗是如何做到比矩形窗更高的副瓣衰减的？

(2)窗函数序数序列的长度N是否会影响窗函数的副瓣衰减?试从理论上加以证明。

(3)通过实验分析并说明凯瑟窗序列中参数β参数与窗函数频率特性之间的关系。

(4）在对信号进行频域分析前通常需要进行加窗处理。截取要求的点数,即相关于加矩窗。也可将被处理信号序列与其他各种类型的窗函数进行相乘后再进行频域分析。通过理论分析或者仿真实例说明加窗处理对谱分析的频谱泄漏效应和分辨力特性有何影响?

(5)常用的矩形窗、汉宁窗、海明窗、布莱克曼窗和凯瑟窗的时域图形是否是中心对称的,这几个窗函数的相频特性是否为线性相位的?混合窗函数的时域图形是否中心对称的，其相频特性是否为线性相位的？

六、总结

通过“窗函数法设计FIR数字滤波器”的实验，我对于滤波器的设计又了更加深入的理解。特别是对于滤波器类型的选择，还有根据所给出的条件进行滤波器的确定还有滤波器窗口长度的确定。对于不同种类的窗口，如：矩形窗、三角窗、汉宁窗、海明窗等都更加熟悉，也知道了它们的表达式和性能特点，以及它们的函数图形。我也通过实验，知道了已知单位冲激响应，要在Matlab中求其幅频特性，并进行作图的方法。对于FIR滤波器的知识也有了更扎实的掌握。对我理解窗函数设计滤波器方法很有帮助，我相信，这次数字信号处理实验的参与仅仅是打开了我利用matlab信号分析与处理的大门，在今后我一定会不断努力加油让matlab真正成为我的一项技能。