【数字信号处理2】IIR 滤波器设计

一、实验目的

1.掌握冲激响应法和双线性变换法设计IIR滤波器的原理及具体设计方法，熟悉用双线性设计法设计低通、带通和高通IIR数字滤波器的计算机程序；
2.熟悉模拟Butterworth滤波器的设计，掌握冲激响应法和双线性变换法设计数字IIR滤波器的方法。

二、实验内容

1、不同阶次模拟巴特沃兹滤波器的频率响应

结论：不同阶次的所对应的滤波器的幅度大值走向一样，但其过渡带存在明显的不同，阶次越高，滤波器的过渡带越小。
2、根据模拟滤波器指标，设计低通巴特沃兹滤波器
设计指标为：通带截止频率fp=6kHz,通带最大衰减ɑp =3dB，阻带截止频率fs= 14kHz,阻带最小衰减ɑs=32dB。
归一化之前的低通巴特沃兹滤波器的频率响应:

归一化之后的低通巴特沃兹滤波器德品频率响应:

3、已知模拟滤波器传输函数H(s)为

按照不同采样频率Fs1=1.1Hz, Fs2=10.1Hz 使用脉冲响应不变法将H(s)转换为H(z)数字IIR滤波器。
不同采样采样频率对应的IIR滤波器的频率响应见下：
结论：不同的采样频率设计出的IIR滤波器的频率响应是明显不一样的，采样频率越高，所设计出的频率计在pi处衰减的越大，再滤波器越接近模拟滤波器。

4 给定数字IIR指标，用脉冲响应不变法设计IIR。
在不同的Fs下，所得到的模拟滤波器和数字滤波器的频率响应见下：

结论：在不同的采样频率下，有数字指标所得到模拟指标不一样，故导致所设计的模拟滤波器的频率响应不一样，但是，最后由模拟滤波器转化为的数字滤波器差别不大，是符合所给的数字指标的。

5、用双线性变换法设计IIR，并同脉冲响应不变法比较设计的结果。
数字滤波器指标要求如下：

将所得的数字滤波器的频率响应与题目给的要求对应，可以得之所设计的滤波器符合题目的要求。

三、实验结论

冲激响应不变法设计的滤波器在直流附近的数字域与模拟域的幅频特性相近，但双线性变换法是以严重的失真换取了频率的不混叠，这一般用于具有片段常数特性的滤波器。

附-实验源码

%1 利用Matlab系统函数计算不同阶次的巴特沃兹滤波器频率响应
clc,clear all,close all;
freq_axis =[0:0.01:2];%归一化的频率
N_array =[461113];%滤波器的阶次
figure(1);for I=1:4
    N =N_array(I)[z p k]=buttap(N)%获得6阶巴特沃兹滤波器的参数
    [b a]=zp2tf( z, p, k )%得到传递函数 tf:trans-function
    disp(['N='num2str(N)'极点:']);
    p
    disp(['N='num2str(N),'零点:']);
    z
    
    [H ]=freqs(b, a, freq_axis);
    magh =abs(H);subplot(2,2, I );plot( freq_axis, magh );axis([0201]);%限制坐标轴
    xlabel('w/wc -- 归一化频率');ylabel('幅度谱');title(['巴特沃兹滤波器幅度谱 N=',num2str(N)]);
    
end

%2 根据提出模拟IIR的指标，计算IIR的阶次、系数等

J =sqrt(-1);
fp =6e3;%6KHz
alfa_p =3;%最大通带衰减3dB
fs =14e3;%14KHz 阻带截止频率
alfa_s =32;%32dB

k_sp =sqrt(10^(alfa_p/10)-1)/sqrt(10^(alfa_s/10)-1)
lamda_sp =2*pi*fs/(2*pi*fp)
N =ceil(-log(k_sp)/log(lamda_sp))

p =zeros(1,N);for I=1:N
    k = I-1;p(I)=exp( J*pi/2+ J*pi*(2*k+1)/2/N );
end
z =[];%表示无零点
[b1 a1]=zp2tf( z, p,1);%得到传递函数 tf:trans-function

fc_freq =  fp*(10^(alfa_p/10)-1)^(-1/2/N);%3dB截至频率
figure(2);
freq_axis1 =[0:50:15e3];[H1 ]=freqs(b1, a1, freq_axis1/fc_freq);%对Wc作频率归一化
plot( freq_axis1 ,20*log10(abs(H1)));

H1 =freqs(b1,a1,freq_axis);figure(3);plot( freq_axis,abs(H1));
grid on

% 给出的模拟滤波器，使用脉冲响应不变法来设计数字IIR滤波器
% 验证不同采样频率fs对于数字IIR滤波器的影响
clc,clear all,close all;
J =sqrt(-1);
b =[000.5250];%高次项系数在前
a =[10.55500.7552];
freq_analog =[0:0.5:30];%模拟角频率!!! 
H_analog =freqs(b, a, freq_analog);%计算在模拟角频率处的频响
figure(1);subplot(2,1,1);plot( freq_analog,20*log10(abs(H_analog)));xlabel('角频率( rad/s )');ylabel('幅频响应(dB)');title('模拟滤波器频率响应');

Fs1 =1.1;%1.1Hz采样率
[bz1, az1]=impinvar( b, a, Fs1 );%根据模拟滤波器系数b,a，按照采样率fs转换为系数为bz,az的数字IIR滤波器

Fs2 =10.1;%10.1Hz采样频率
[bz2, az2]=impinvar( b, a, Fs2 );%根据模拟滤波器系数b,a，按照采样率fs转换为系数为bz,az的数字IIR滤波器

freq_digital =[0:0.1:2*pi];
H_dig1 =freqz( bz1, az1, freq_digital );
H_dig2 =freqz( bz2, az2, freq_digital );subplot(2,1,2);plot( freq_digital,20*log10(abs(H_dig1)));xlabel('数字频率');ylabel('幅频响应(dB)');title('不同采样频率对于数字IIR频响的影响');
hold on;plot( freq_digital,20*log10(abs(H_dig2)),'r');legend(['Fs1='num2str(Fs1)'Hz'],['Fs2='num2str(Fs2)'Hz']);

clc,clear all,close all;%1 给定数字滤波器指标，用脉冲响应不变法设计数字滤波器
J =sqrt(-1);
freq_p =0.2*pi;%通带截止频率0.2pi
alpha_p =1;%通带最大衰减1dB
freq_s =0.3*pi;
alpha_s =15;

Fs_array =[110];
Color_array =['b','k'];for III=1:2
    color_used =Color_array(III);
    Fs =Fs_array(III);
    T =1/Fs;
    analog_freq_p = freq_p/T   %数字频率w到模拟角频率W的转换 W = w/T;
    analog_freq_s = freq_s/T

    %根据模拟滤波器的指标求巴特沃兹滤波器的阶次
    k_sp =sqrt(10^(alpha_p/10)-1)/sqrt(10^(alpha_s/10)-1);
    lamda_sp =2*pi*analog_freq_s/(2*pi*analog_freq_p);
    N =ceil(-log(k_sp)/log(lamda_sp))
    Wc =  analog_freq_p*(10^(analog_freq_p/10)-1)^(-1/2/N)%3dB截至频率
   

    %查表，得到归一化的传输函数 或者根据公式得到极点p
    p =zeros(1,N);for I=1:N
        k = I-1;p(I)=exp( J*pi/2+ J*pi*(2*k+1)/2/N );
    end

    % 根据得到的零点、极点得到归一化的传输函数H(p)
    z =[];%表示无零点
    [b a]=zp2tf( z, p,1);%得到传递函数 tf:trans-function,第三个参数为增益=1% 根据p=s/Wc 把归一化传输函数H(p)的系数转换为实际的模拟传输函数H(s)
    bs =zeros(1,length(b))%H(s)的分子系数
    as =zeros(1,length(a))% 对与N阶模拟滤波器，设其传输函数分子分母系数分别为bs,as
    temp =[N:-1:0]
    temp = Wc .^ temp
    bs = b ./ temp
    as = a ./ temp

    temp =as(1)
    as = as/temp             %验证与书上结果是否一致 
    bs = bs/temp
    figure(1);
    freq_analog =linspace(0,2*analog_freq_s );
    H_analog =freqs(bs,as,freq_analog);plot( freq_analog,20*log10(abs(H_analog)),color_used); grid on; hold on;xlabel('模拟角频率(rad/s)');ylabel('幅频响应(dB)');title('模拟滤波器的频率响应');% 根据H(s)获得H（z），使用脉冲响应不变法
    [bz, az]=impinvar( bs, as, Fs );%根据模拟滤波器系数b,a，按照采样率fs转换为系数为bz,az的数字IIR滤波器
    freq_digital =[0:0.1:2*pi-0.1];
    H =freqz( bz, az, freq_digital );figure(2);plot( freq_digital/pi,20*log10(abs(H)),color_used ); grid on; hold on;xlabel('数字频率/pi');ylabel('幅频响应(dB)');title('数字IIR滤波器频率响应');
    
    a =1;

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 当采样频率变为2Hz时，求得出的数字IIR频率响应
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%figure(1),legend(['Fs='num2str(Fs_array(1))'Hz'],['Fs='num2str(Fs_array(2))'Hz']);figure(2),legend(['Fs='num2str(Fs_array(1))'Hz'],['Fs='num2str(Fs_array(2))'Hz']);

clc,clear all,close all;%  给定数字滤波器指标
J =sqrt(-1);
freq_p =0.3*pi;%通带截止频率0.3pi
alpha_p =1.1;%通带最大衰减1.1dB
freq_s =0.5*pi;
alpha_s =16;% 模拟低通指标
Fs =1;
T =1/Fs;
analog_freq_p =2/T *tan(freq_p/2)
analog_freq_s =2/T *tan(freq_s/2)% 计算巴特沃兹滤波器的指标
k_sp =sqrt(10^(alpha_p/10)-1)/sqrt(10^(alpha_s/10)-1);
lamda_sp =2*pi*analog_freq_s/(2*pi*analog_freq_p);
N =ceil(-log(k_sp)/log(lamda_sp));
Wc =  analog_freq_p*(10^(analog_freq_p/10)-1)^(-1/2/N);%3dB截至频率
   
%查表，得到归一化的传输函数 或者根据公式得到极点p
p =zeros(1,N);for I=1:N
    k = I-1;p(I)=exp( J*pi/2+ J*pi*(2*k+1)/2/N );
end

% 根据得到的零点、极点得到归一化的传输函数H(p)
z =[];%表示无零点
[b a]=zp2tf( z, p,1);%得到传递函数 tf:trans-function,第三个参数为增益=1% 根据p=s/Wc 把归一化传输函数H(p)的系数转换为实际的模拟传输函数H(s)
bs =zeros(1,length(b));%H(s)的分子系数
as =zeros(1,length(a));% 对与N阶模拟滤波器，设其传输函数分子分母系数分别为bs,as
temp =[N:-1:0];
temp = Wc .^ temp;
bs = b ./ temp;
as = a ./ temp;

temp =as(1);
as = as/temp;%验证与书上结果是否一致 
bs = bs/temp;figure(1);
freq_analog =linspace(0,2*analog_freq_s );
H_analog =freqs(bs,as,freq_analog);plot( freq_analog,20*log10(abs(H_analog)),'k'); grid on;xlabel('模拟角频率(rad/s)');ylabel('幅频响应(dB)');title('模拟滤波器的频率响应');% 根据模拟滤波器的传输函数得到数字滤波器的传输函数
[bz az]=bilinear( bs, as, Fs );% 当模拟滤波器bs,as为其系数，用Fs采样。转换成数字IIR的传输函数bz,az
freq_digital =[0:0.1:2*pi-0.1];
H =freqz( bz, az, freq_digital );figure(2);plot( freq_digital/pi,20*log10(abs(H)),'k'); grid on;xlabel('数字频率/pi');ylabel('幅频响应(dB)');title('数字IIR滤波器频率响应');% 可以比较一下用脉冲响应不变法获得频率响应