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 = [ 1 2 ];
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'] );
