hc822.m

这是数字信号处理——matlab释义与实现一书的源代码

%《数字信号处理教程——MATLAB释义与实现》第八章例8.2.2程序hc822
% 巴特沃斯模拟滤波器设计
% 电子工业出版社出版  陈怀琛编著 2004年9月
%
OmegaP=2*pi*5000;OmegaS=2*pi*12000; % 给出原始要求
Rp=1;As=30;
N=ceil( log10((10.^(0.1*abs(Rp))-1)./(10.^(0.1*abs(As))-1))/(2*log10(OmegaP/OmegaS)));  % 估算N的设计公式
OmegaC= OmegaP / ( (10^(.1*abs(Rp)) - 1)^(1/(2*N)) );   % 估算OmegaC的设计公式
%[N, OmegaC] =buttord(OmegaP, OmegaS, Rp, As, 's')      % 调用信号处理工具箱中估算 N 和 OmegaC 的公式
[z0,p0,k0] = buttap(N)          % 调用buttap函数,得到系统的左半平面零极点
p = p0*OmegaC; z = z0*OmegaC;   		% 将零极点乘以Omegac，得到非归一化零极点
b0=k0*real(poly(z0)), a0=real(poly(p0))         % 求归一化滤波器系数b0,a0   
Nm=length(a0)-length(b0);		% 求极点个数与零点个数之差Nm
k = k0*OmegaC^Nm;			 		% 将k0乘以Omegac^Nm,得到非归一化k
b=k*real(poly(z)); a=real(poly(p))         % 求非归一化滤波器系数b,a   
pause
wx=[OmegaP,OmegaS]/OmegaC;              % 设置频率向量[OmegaP,OmegaS]
Hx=freqs(b0,a0,wx)                      % 计算该两点上的幅特性
dbHx=20*log10(abs(Hx))                  % 化为分贝值
freqs(b,a)                              % 画出波特图
[G,sos]=tf2sos(b0,a0)                   % 化为二阶级联形式
set(gcf,'color','w')                  % 置图形背景色为白