河海大学通信工程数字信号实验5.docx
《河海大学通信工程数字信号实验5.docx》由会员分享，可在线阅读，更多相关《河海大学通信工程数字信号实验5.docx（22页珍藏版）》请在冰点文库上搜索。
河海大学通信工程数字信号实验5
实验五有限长单位脉冲响应滤波器设计
一、实验目的
（1）掌握用窗函数法、频率采样法及优化设计法设计FIR滤波器的原理及方法，熟悉相应的MATLAB编程。
（2）熟悉线性相位FIR滤波器的幅频特性和相频特性。
（3）了解各种不同窗函数对滤波器性能的影响。
二、实验原理
window=ones（1,N）:
产生N点矩形窗，行向量。
window=hann（N）:
产生N点汉宁窗，列向量。
window=hanning（N）:
产生N点非零汉宁窗，列向量。
等价于去除hann（N+2）的第一个零元素和最后一个零元素，得到的N点非零窗函数。
window=hamming（N）:
产生N点海明窗，列向量。
window=blackman（N）:
产生N点布莱克曼窗，列向量。
window=kaiser（N,beta）:
产生参数为beta的N点凯塞窗，列向量。
[M,Wd,beta,ftype]=kaiserord（f,a,dev,fs）:
凯塞窗参数估计。
f为一组边界频率，最高频率为fs/2。
a为f中各个频带的幅度值，通带取1，阻带取0。
如果f中有2个元素，则形成3个频带，其中第1个和第3个是通带或阻带，第2个是过渡带，a中也有2个元素，指明第1个和第3个频带是通带还是阻带；如果f中有4个元素，则形成5个频带，其中1，3和5是通带或阻带，2和4是过渡带，a中有3个元素，指明1，3和5是通带还是阻带。
dev的维数与a相同，指明每个频带上的波动值。
fs为采样频率。
M为FIR滤波器的阶数，M=N-1。
Wd为归一化边界频率，等于数字边界角频率除以π，或者边界频率除以fs/2。
beta就是凯塞窗的参数β。
ftype为滤波器的类型。
fs=8000;
f=[10001500];
a=[10];
dev=[0.050.01];
[M,Wd,beta,ftype]=kaiserord（f,a,dev,fs）;
h=fir1（M,Wd,ftype,kaiser（M+1,beta））;
freqz（h）;
fs=8000;
f=[1000130022102410];
a=[010];
dev=[0.010.050.01];
[M,Wd,beta,ftype]=kaiserord（f,a,dev,fs）;
hh=fir1（M,Wd,ftype,kaiser（M+1,beta））;
[H,f]=freqz（hh,1,1024,fs）;
figure;plot（f,abs（H））;gridon;
b=fir1（M,Wd,'ftype',window）:
用窗函数法求FIR滤波器的系数b（单位脉冲响应）。
M为滤波器的阶数，M=N-1。
Wd为一组归一化边界频率，通带和阻带间隔分布，无过渡带；只有一个元素，表示低通或高通滤波器；有两个元素表示带通和带阻滤波器；有三个及以上元素，表示多带滤波器。
'ftype'表示滤波器类型，'high'表示高通滤波器，'stop'表示带阻滤波器，'DC-0'表示多带滤波器的第一个频带为阻带，'DC-1'表示多带滤波器的第一个频带为通带。
window为窗口类型，缺省为海明窗。
b=fir2（M,f,m,window）:
用频率采样法求FIR滤波器的系数b。
M为滤波器的阶数，M=N-1。
f为一组归一化频率，第一个元素必须为0，最后一个元素必须为1（对应奈奎斯特频率，即采样频率的一半），中间的元素按升序排列。
m的维数与f相同，指明f中每个频率上的理想幅度。
window为窗口类型，缺省为海明窗。
Fir2可以实现任意幅度特性的滤波器。
三、实验内容
1、用窗函数法设计一个线性相位FIR高通滤波器，通带边界频率为0.6π，阻带边界频率为0.4π，要求阻带衰减不小于40dB，窗函数从矩形窗、汉宁窗、海明窗和布莱克曼窗中选取，且要求滤波器的阶数最小。
给出该滤波器的单位脉冲响应h（n）的解析式，并用MATLAB绘出时域波形和幅频特性。
2、用窗函数法设计一个线性相位FIR带通滤波器，通带边界频率为4kHz和6kHz，阻带边界频率为2kHz和8kHz，采样频率为20kHz，要求阻带衰减不小于50dB，窗函数从矩形窗、汉宁窗、海明窗和布莱克曼窗中选取，且要求滤波器的阶数最小。
给出该滤波器的单位脉冲响应h（n）的解析式，并用MATLAB绘出时域波形和幅频特性。
3、用窗函数法设计一线性相位FIR带阻滤波器，通带边界频率为600Hz和1400Hz，阻带边界频率为800Hz和1200Hz，采样频率为4000Hz，要求阻带衰减不小于50dB，窗函数从矩形窗、汉宁窗、海明窗和布莱克曼窗中选取，且要求滤波器的阶数最小。
给出该滤波器的单位脉冲响应h（n）的解析式，并用MATLAB绘出时域波形和幅频特性。
4、用凯塞窗设计一个多带线性相位滤波器，幅频特性如下图所示，N=40，β分别取4、6、10，比较不同β值时的幅频特性和相频特性。
5、用频率采样法设计第（4）题的多带滤波器，比较过渡带不设采样点和设置一个采样点的幅频特性，过渡点的值设为0.5，窗函数选矩形窗。
四、实验结果与分析
1、用窗函数法设计一个线性相位FIR高通滤波器，通带边界频率为0.6π，阻带边界频率为0.4π，要求阻带衰减不小于40dB，窗函数从矩形窗、汉宁窗、海明窗和布莱克曼窗中选取，且要求滤波器的阶数最小。
给出该滤波器的单位脉冲响应h（n）的解析式，并用MATLAB绘出时域波形和幅频特性。
解：
由题可得:
求理想低通滤波器的边界频率。
理想低通滤波器的边界频率是实际滤波器幅度函数正负肩峰频率的中心，但这里不知道正负肩峰的频率，因此用通带和阻带边界频率的中心近似:
选择窗函数和窗口长度。
阻带衰减不小于40dB，因此选择海明窗（其阻带最小衰减为-44dB）；滤波器的过渡带宽为0.6π-0.4π=0.2π，因此窗口长度为:
线性相位延迟常数为:
根据理想边界频率
和线性相位延迟常数
，求理想单位脉冲响应
:
窗函数与理想单位脉冲响应相乘，得到线性相位FIR低通滤波器的单位脉冲响应:
matlab程序见附录程序一
运行结果如下：
分析：
由图知阻带衰减最小值大于40，满足要求。
2、用窗函数法设计一个线性相位FIR带通滤波器，通带边界频率为4kHz和6kHz，阻带边界频率为2kHz和8kHz，采样频率为20kHz，要求阻带衰减不小于50dB，窗函数从矩形窗、汉宁窗、海明窗和布莱克曼窗中选取，且要求滤波器的阶数最小。
给出该滤波器的单位脉冲响应h（n）的解析式，并用MATLAB绘出时域波形和幅频特性。
解：
求数字边界频率:
求理想滤波器的边界频率:
求理想单位脉冲响应:
选择窗函数。
阻带最小衰减为-50dB，因此选择汉明窗（其阻带最小衰减为-53dB）；滤波器的过渡带宽为0.4π-0.2π=0.8π-0.6π=0.2π，因此
求FIR滤波器的单位脉冲响应h（n）:
matlab程序见附录程序二
运行结果如下：
分析：
由图可知，阻带衰减不小于50dB满足要求。
3、用窗函数法设计一线性相位FIR带阻滤波器，通带边界频率为600Hz和1400Hz，阻带边界频率为800Hz和1200Hz，采样频率为4000Hz，要求阻带衰减不小于50dB，窗函数从矩形窗、汉宁窗、海明窗和布莱克曼窗中选取，且要求滤波器的阶数最小。
给出该滤波器的单位脉冲响应h（n）的解析式，并用MATLAB绘出时域波形和幅频特性。
解：
求数字边界频率:
求理想滤波器的边界频率:
求理想单位脉冲响应:
（1）选择窗函数。
阻带最小衰减为-50dB，因此选择汉明窗（其阻带最小衰减为-53dB）；滤波器的过渡带宽为0.4π-0.3π=0.7π-0.6π=0.1π，因此
求FIR滤波器的单位脉冲响应h（n）:
matlab程序见附录程序三
运行结果如下：
分析：
由图可知，阻带衰减不小于50dB，满足要求。
4、用凯塞窗设计一个多带线性相位滤波器，幅频特性如下图所示，N=40，β分别取4、6、10，比较不同β值时的幅频特性和相频特性。
解：
matlab程序见附录程序四
运行结果如下：
分析：
β=4时幅频图比较接近理想滤波器，相频特性不太明显，随着β的增大，幅频图过渡带开始增加，相频特性越来越明显。
这可能是因为增大β导致频谱的主瓣宽度增加造成的。
5、用频率采样法设计第（4）题的多带滤波器，比较过渡带不设采样点和设置一个采样点的幅频特性，过渡点的值设为0.5，窗函数选矩形窗。
解：
由于（2π/40）×4=0.2π，（2π/32）×8=0.4π，（2π/40）×12=0.6π，（2π/40）×16=0.8π，又因为N为偶数，属于第二类线性相位滤波器，幅度函数在0~2π上以π为中心呈奇数对成分布。
因此频率采样值的幅度为:
频响采样的相位为:
matlab程序见附录程序五
运行结果如下：
分析：
在理想频率响应的间断点的边缘加上一个过渡采样点，阻带衰减变大，但过渡带也就越宽。
五、实验小结
通过此次实验，掌握了用窗函数法、频率采样法及优化设计法设计FIR滤波器的原理及方法，熟悉了相应的MATLAB编程。
熟悉了线性相位FIR滤波器的幅频特性和相频特性。
了解了各种不同窗函数对滤波器性能的影响。
在实验中，明白了matlab在数字信号处理中的重要性，学好matlab对以后工作学习有极大的帮助。
六、附录
程序一
clear;
N=31;n=0:
N-1;
hd=（sin（pi*（n-15））-sin（0.5*pi*（n-15）））./（pi*（n-15））;hd（16）=0.5;
win=hanning（N）;
h=win'.*hd;
figure;stem（n,h）;
xlabel（'n'）;ylabel（'h（n）'）;grid;
title（'FIR高通滤波单位脉冲响应h（n）'）;
[H,w]=freqz（h,1）;H=20*log10（abs（H））;
figure;
plot（w/pi,H）;
axis（[01-10010]）;
xlabel（'\omega/\pi'）;ylabel（'幅度/dB'）;grid;
title（'FIR高通滤波器，hanning窗，N=31'）;
程序二
clear;
N=33;n=0:
N-1;
hd=（sin（0.7*pi*（n-16））-sin（0.3*pi*（n-16）））./（pi*（n-16））;hd（17）=0.4;
win=hamming（N）;
h=win'.*hd;
figure;stem（n,h）;
xlabel（'n'）;ylabel（'h（n）'）;grid;
title（'FIR带通滤波单位脉冲响应h（n）'）;
[H,w]=freqz（h,1）;H=20*log10（abs（H））;
figure;
plot（w/pi,H）;
axis（[01-10010]）;
xlabel（'\omega/\pi'）;ylabel（'幅度/dB'）;grid;
title（'FIR带通滤波器，hamming窗，N=33'）;
程序三
clear;
N=66;n=0:
N-1;
hd=（sin（pi*（n-32.5））+sin（0.35*pi*（n-32.5））-sin（0.65*pi*（n-32.5）））./（pi*（n-32.5））;
win=hamming（N）;
h=win'.*hd;
figure;stem（n,h）;
xlabel（'n'）;ylabel（'h（n）'）;grid;
title（'FIR带通滤波单位脉冲响应h（n）'）;
[H,w]=freqz（h,1）;H=20*log10（abs（H））;
figure;
plot（w/pi,H）;
axis（[01-10010]）;
xlabel（'\omega/\pi'）;ylabel（'幅度/dB'）;grid;
title（'FIR带阻滤波器，hamming窗，N=66'）;
程序四
clear;
%β=4
Wd=[0.20.40.60.8];M=39;beta=4;
hh=fir1（M,Wd,'DC-0',kaiser（M+1,beta））;
[H,w]=freqz（hh,1）;
figure;subplot（2,1,1）;
plot（w/pi,abs（H））;
xlabel（'\omega/\pi'）;ylabel（'幅度'）;grid;
title（'幅频特性,beta=4'）;
subplot（2,1,2）;
plot（w/pi,angle（H））;
xlabel（'\omega/\pi'）;ylabel（'弧度'）;grid;
title（'相频特性'）;
%β=6;
Wd=[0.20.40.60.8];M=39;beta=6;
hh=fir1（M,Wd,'DC-0',kaiser（M+1,beta））;
[H,w]=freqz（hh,1）;
figure;subplot（2,1,1）;
plot（w/pi,abs（H））;
xlabel（'\omega/\pi'）;
ylabel（'幅度'）;grid;
title（'幅频特性,beta=6'）;
subplot（2,1,2）;
plot（w/pi,angle（H））;
xlabel（'\omega/\pi'）;
ylabel（'弧度'）;grid;
title（'相频特性'）;
%β=10;
Wd=[0.20.40.60.8];M=39;beta=10;
hh=fir1（M,Wd,'DC-0',kaiser（M+1,beta））;
[H,w]=freqz（hh,1）;
figure;subplot（2,1,1）;
plot（w/pi,abs（H））;
xlabel（'\omega/\pi'）;ylabel（'幅度'）;grid;
title（'幅频特性,beta=10'）;
subplot（2,1,2）;
plot（w/pi,angle（H））;
xlabel（'\omega/\pi'）;ylabel（'弧度'）;grid;
title（'相频特性'）;
程序五
clear;
N=40;
k=0:
N-1;
%过渡带不设采样点
Hk=[zeros（1,4）,ones（1,5）,zeros（1,3）,ones（1,5）,zeros（1,7）,-ones（1,5）,zeros（1,3）,-ones（1,5）zeros（1,3）];
h=real（ifft（Hk.*exp（-j*pi*（N-1）*k/N）））;
[H,w]=freqz（h,1）;
figure;plot（w/pi,20*log10（abs（H）））;
axis（[01-7010]）;grid;
xlabel（'\omega/\pi'）;ylabel（'幅度/dB'）;
title（'频率采样法，N=40，过渡带不设采样点'）;
%过渡带设一个采样点
Hk=[zeros（1,3）0.5ones（1,5）0.5zeros（1,1）0.5ones（1,5）0.5zeros（1,5）-0.5-ones（1,5）-0.5zeros（1,1）-0.5-ones（1,5）-0.5zeros（1,2）];
h=real（ifft（Hk.*exp（-j*pi*（N-1）*k/N）））;
[H,w]=freqz（h,1）;
figure;plot（w/pi,20*log10（abs（H）））;
axis（[01-7010]）;grid;
xlabel（'\omega/\pi'）;ylabel（'幅度/dB'）;
title（'频率采样法，N=40，过渡带设一个采样点'）;
升级会员 