hc835.m

数字信号处理的matlab程序

%《数字信号处理教程——MATLAB释义与实现》第八章例8.3.5程序hc835
% 双线性变换法由巴特沃思变数字滤波器
% 电子工业出版社出版  陈怀琛编著 2004年9月
%
% 数字滤波器指标:
wp = 0.2*pi;                         % 数字通带频率(弧度)
ws = 0.3*pi;                         % 数字阻带频率(弧度)
Rp = 1;                              % 通带波动(dB)
As = 15;                             % 阻带衰减(dB) 

% 模拟原型指标的频率逆映射
T = 0.001; Fs = 1/T;                     % 置 T=0.001
OmegaP = (2/T)*tan(wp/2);            % 原型通带频率预修正
OmegaS = (2/T)*tan(ws/2);            % 原型阻带频率预修正
ep = sqrt(10^(Rp/10)-1);             % 通带波动参数
Ripple = sqrt(1/(1+ep*ep));          % 通带波动
Attn = 1/(10^(As/20));               % 阻带衰减

% 模拟巴特沃思原型滤波器计算:
[N,OmegaC] =buttord(OmegaP,OmegaS,Rp,As,'s');   % 原型的阶数和截止频率计算
%%***巴特沃思滤波器阶次 =  6 
[z0,p0,k0] = buttap(N);         % 归一化巴特沃思原型设计函数
p = p0*OmegaC; z = z0*OmegaC;   % 将零极点乘以Omegac，得到非归一化零极点
k = k0*OmegaC^N;                % 将k0乘以Omegac^N,得到非归一化k
ba0 = real(poly(z0));ba0 = k0*ba0       % 由零点计算分子系数向量
aa0 = real(poly(p0))               % 由极点计算分母系数向量
ba = real(poly(z));ba = k*ba       % 由零点计算分子系数向量
aa = real(poly(p))               % 由极点计算分母系数向量
[bd,ad] = bilinear(ba,aa,Fs); 	% 双线性变换:
[bd1,ad1] = bilinear(ba0,aa0,Fs/OmegaC); 	% 双线性变换:
[sos,G] = tf2sos(bd,ad)			    % 变为二阶环节级联结构    

% 绘图
figure(1); subplot(1,1,1)
[db,mag,pha,grd,w] = myfreqz(bd1,ad1);               % 检验频率响应
subplot(2,2,1); plot(w/pi,mag); title('幅度响应')
xlabel(''); ylabel('|H|');axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]); grid    % 画刻度线
%set(gca,'YTickmode','manual','YTick',[0,Attn,Ripple,1])
subplot(2,2,3); plot(w/pi,db); title('模值(dB)');
xlabel('频率：（单位：pi）'); ylabel('分贝'); axis([0,1,-40,5]);
set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]);        % 画刻度线
set(gca,'YTickmode','manual','YTick',[-50,-15,-1,0]); grid
set(gca,'YTickLabelMode','manual','YTickLabels',['50';'15';' 1';' 0'])
subplot(2,2,2); plot(w/pi,pha/pi); title('相位响应')
xlabel(''); ylabel('单位:pi'); axis([0,1,-1,1]);      
set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]);        % 画刻度线
set(gca,'YTickmode','manual','YTick',[-1,0,1]); grid
subplot(2,2,4); plot(w/pi,grd); title('群延迟')
xlabel('频率（单位：pi）'); ylabel('样本'); axis([0,1,0,10])
set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]);        % 画刻度线
set(gca,'YTickmode','manual','YTick',[0:2:10]); grid
set(gcf,'color','w')              % 置图形背景色为白