【数字信号处理3】FIR滤波器设计

一、实验目的

1.熟悉FIR滤波器设计的基本方法；
2.掌握窗函数法设计FIR滤波器的原理和方法；
3.熟悉线性相位FIR滤波器的幅频特性和相位特性；
4.了解不同窗函数对滤波器性能的响应。

二、实验内容

1、验证窗函数N变化时，验证其频谱主瓣副瓣比、主瓣宽度的变化。
矩形窗函数的N变化时，验证其其频谱主瓣副瓣幅度比基本不会发生变化，而主瓣宽度将会变窄。这说明，当用矩形窗函数设计滤波器时，增大N不能使得阻带衰减减小，但能够减小过渡带。
取N分别为24和34，其频率响应见下：
N为24 时：

由图可以算得：频谱主瓣副瓣幅度比为1.3409，过渡带为：0.48pi，约等于4*pi/N。
N为34 时：

由图可以算得：频谱主瓣副瓣幅度比为1.3719，过渡带为：0.39pi，约等于4*pi/N。
结论：矩形窗函数的N变化时，验证其其频谱主瓣副瓣幅度比基本不会发生变化，而主瓣宽度将会变窄。这说明，当用矩形窗函数设计滤波器时，增大N不能使得阻带衰减减小，但能够减小过渡带。
2、用窗函数法设计线性相位FIR低通
通带截止频率wp=0.49PI, 阻带截止频率ws=0.59PI, 阻带衰减不小于40dB,通带衰减不大于3dB.
选取Hanning和Hamming就可以满足设计的要求。
根据过渡带的要求选择窗的长度。过渡带是0.1 PI,海明和汉宁窗的主瓣宽度均是8 PI/N的，从而选择的窗的长度都是大于等于80的，选择N大于等于80就可以了，在这里选取N=81.
用Hanning和Hamming窗设计的滤波器见下面：

从图中可以看出，这两者都符合设计要求的指标。
3、用窗函数法设计线性相位FIR高通，通带截止频率wp=0.79PI, 阻带截止频率ws=0.69PI, 阻带衰减不小于30dB,通带衰减不大于3dB.
1）根据阻带衰减的要求选择窗函数，可以选择汉宁窗（最大衰减是44dB）
2）根据过渡带的要求选择窗的长度。过渡带是0.1 PI,海明和汉宁窗的主瓣宽度均是8 PI/N的，从而选择的窗的长度都是大于等于80的，选择N=81（必须是奇数）;

由图可以看出，此设计符合指标要求。
结论：窗的长度必须是奇数，而不能是偶数。

三、实验结论

FIR滤波器的特点是阶数较高，但是只要h（n）是实序列，且关于(N-1)/2对称就可以保证严格的线性相位。
一般来说采样窗函数特点是简单，容易实现。但是窗函数不容易进行边界频率的控制，而且改变参数对于幅频特性的影响是不够直观。可以采用频率抽样结构来实现。而且可知频率抽取结构的特点是利用零极点相互抵消来实现的。这样就容易造成系统的不稳定。一般可以采用对系统的修正使得极点都在单位圆内的，而且要采用过渡带内插法来增大阻带的衰减。

附-实验代码

 N =24;
n1=ones(1,N);
fft_1=fft(n1,128);
fft_1=fft_1/max(fft_1);
 
x_axis=(0:1:127)/128*2*3.14;
figure (1);
plot(x_axis,20*log10(abs(fft_1)));grid on;
title('N=24')


 N =34;
n1=ones(1,N);
fft_1=fft(n1,128);
fft_1=fft_1/max(fft_1);
x_axis=(0:1:127)/128*2*3.14;
figure(2) 
plot(x_axis,20*log10(abs(fft_1)));grid on;
title('N=34')

clc,clear all,close all;
wp = 0.49*pi;        %通带截至频率
ws = 0.59*pi;        %阻带截至频率
wdel = ws - wp;     %过渡带宽

% Hanning窗--------------------------------------------------------------
N_hanning = ceil( 8*pi/wdel );
Wn = (wp + ws)/2   %截止频率

% N_hanning取奇数
if mod(N_hanning,2)==0
    N_hanning = N_hanning + 1;
end
window_hanning = hanning(N_hanning)   %获得hanning窗离散序列

b_hanning = fir1( N_hanning-1, Wn/pi, window_hanning )  
%指定滤波器阶次，归一化截止频率，窗函数序列得到h[n]
%注意滤波器的阶次=窗长-1   !!!!!

freq_axis = [0:pi/512:pi-pi/512];   
freq_norm = [0:511]/512;    %归一化的频率轴     
H_hanning = freqz( b_hanning, 1, 512);
subplot(2,1,1);
plot( freq_norm,20*log10(abs(H_hanning)) ); hold on;
xlabel( '归一化频率w/pi' ); ylabel( '幅度(dB)' );
title( '采用hanning和hamming设计的FIR-幅度响应');

subplot(2,1,2);
plot( freq_norm,angle(H_hanning) ); hold on;
xlabel( '归一化频率w/pi' ); ylabel( '相位' );
title( '采用hanning和hamming设计的FIR-相位响应');



% Hamming窗--------------------------------------------------------------
N_hamming = ceil( 8*pi/wdel );
Wn = (wp + ws)/2;   %截止频率

% N_hamming取奇数
if mod(N_hamming,2)==0
    N_hamming = N_hamming + 1
end
window_hamming = hamming(N_hamming);    %获得hamming窗离散序列

b_hamming = fir1( N_hamming-1, Wn/pi, window_hamming );   %指定滤波器阶次，归一化截止频率，窗函数序列得到h[n]
                                                          %注意滤波器的阶次=窗长-1
H_hamming = freqz( b_hamming, 1, 512);
subplot(2,1,1);
plot(freq_norm, 20*log10(abs(H_hamming)),'k' );

subplot(2,1,2);
plot( freq_norm,angle(H_hamming),'k' ); hold on;


% Hamming窗--------------------------------------------------------------
% 增大N时，查看滤波器带宽以及阻带衰减的变化
N_hamming_2N = ceil( 8*pi/wdel )*2;
Wn = (wp + ws)/2;   %截止频率

% N_hamming_2N取奇数
if mod(N_hamming_2N,2)==0
    N_hamming_2N = N_hamming_2N + 1;
end
window_hamming = hamming(N_hamming_2N);    %获得hamming窗离散序列

b_hamming = fir1( N_hamming_2N-1, Wn/pi, window_hamming );   %指定滤波器阶次，归一化截止频率，窗函数序列得到h[n]
                                                          %注意滤波器的阶次=窗长-1
H_hamming = freqz( b_hamming, 1, 512);
subplot(2,1,1);
plot(freq_norm, 20*log10(abs(H_hamming)),'r' );

subplot(2,1,2);
plot( freq_norm,angle(H_hamming),'r' ); hold on;


legend( ['Hanning' num2str(N_hanning) '阶'] ,['Hamming ' num2str(N_hamming) '阶'],['Hamming ' num2str(N_hamming_2N) '阶'] );

clc,clear all,close all;
wp = 0.79*pi;        %通带截至频率
ws = 0.69*pi;        %阻带截至频率
wdel = wp - ws;     %过渡带宽

% Hanning窗--------------------------------------------------------------
N_hanning = ceil( 8*pi/wdel );
Wn = (wp + ws)/2;   %截止频率

% N_hanning取奇数
if mod(N_hanning,2)==0
    N_hanning = N_hanning + 1;
end
window_hanning = hanning(N_hanning);    %获得hanning窗离散序列

b_hanning = fir1( N_hanning-1, Wn/pi, 'high',window_hanning );  
%指定滤波器阶次，归一化截止频率，窗函数序列得到h[n]
%注意滤波器的阶次=窗长-1   !!!!!

freq_axis = [0:pi/512:pi-pi/512];   
freq_norm = [0:511]/512;    %归一化的频率轴     
H_hanning = freqz( b_hanning, 1, 512);
subplot(2,1,1);
plot( freq_norm,20*log10(abs(H_hanning)) ); hold on;
xlabel( '归一化频率w/pi' ); ylabel( '幅度(dB)' );
title( '采用hanning和hamming设计的FIR-幅度响应');

subplot(2,1,2);
plot( freq_norm,angle(H_hanning) ); hold on;
xlabel( '归一化频率w/pi' ); ylabel( '相位' );
title( '采用hanning和hamming设计的FIR-相位响应');



% Hamming窗--------------------------------------------------------------
N_hamming = ceil( 8*pi/wdel );
Wn = (wp + ws)/2;   %截止频率

% N_hamming取奇数
if mod(N_hamming,2)==0
    N_hamming = N_hamming + 1;
end
window_hamming = hamming(N_hamming);    %获得hamming窗离散序列

b_hamming = fir1( N_hamming-1, Wn/pi, 'high',window_hamming );   %指定滤波器阶次，归一化截止频率，窗函数序列得到h[n]
                                                          %注意滤波器的阶次=窗长-1
H_hamming = freqz( b_hamming, 1, 512);
subplot(2,1,1);
plot(freq_norm, 20*log10(abs(H_hamming)),'k' );

subplot(2,1,2);
plot( freq_norm,angle(H_hamming),'k' ); hold on;


% Hamming窗--------------------------------------------------------------
% 增大N时，查看滤波器带宽以及阻带衰减的变化
N_hamming_2N = ceil( 8*pi/wdel )*2;
Wn = (wp + ws)/2;   %截止频率

% N_hamming_2N取奇数
if mod(N_hamming_2N,2)==0
    N_hamming_2N = N_hamming_2N + 1;
end
window_hamming = hamming(N_hamming_2N);    %获得hamming窗离散序列

b_hamming = fir1( N_hamming_2N-1, Wn/pi,'high', window_hamming );   %指定滤波器阶次，归一化截止频率，窗函数序列得到h[n]
                                                          %注意滤波器的阶次=窗长-1
H_hamming = freqz( b_hamming, 1, 512);
subplot(2,1,1);
plot(freq_norm, 20*log10(abs(H_hamming)),'r' );

subplot(2,1,2);
plot( freq_norm,angle(H_hamming),'r' ); hold on;


legend( ['Hanning' num2str(N_hanning) '阶'] ,['Hamming ' num2str(N_hamming) '阶'],['Hamming ' num2str(N_hamming_2N) '阶'] );