zmapping.m

digital signal processing-MATLAB

function [bz,az] = zmapping(bZ,aZ,Nz,Dz)
%《数字信号处理教程——MATLAB释义与实现》子程序
% 从Z域到z域的频带变换子程序
% 电子工业出版社出版  陈怀琛编著 2004年9月
%
% [bz,az] = zmapping(bZ,aZ,Nz,Dz)
% -------------------------------------------------------
% 实现:
%          b(z)   b(Z)|
%          ---- = ----|     N(z)
%          a(z)   a(Z)|@Z = ----
%                           D(z)
%
%  bz,az = 变换后的滤波器分子分母系数向量
%  bZ,aZ,= 变换前的滤波器分子分母系数向量
%  Nz,Dz = 变换所用的算子的分子分母系数向量
bzord = (length(bZ)-1)*(length(Nz)-1);			% 确定变换后的分子阶数
azord = (length(aZ)-1)*(length(Dz)-1); 		% 确定变换后的分母阶数

bz = zeros(1,bzord+1);				% 分子系数向量初始化为零向量
for k = 0:bzord					% 依次求各多项式乘积结果
    pln = [1];
    for l = 0:k-1  pln = conv(pln,Nz);    end		% 求Nz的k 次乘积
    pld = [1];
    for l = 0:bzord-k-1  pld = conv(pld,Dz);  end	 % 求Dz的bzord-k 次乘积
    bz = bz+bZ(k+1)*conv(pln,pld);			% 求bZ*Nz^k*Dz^(bzord-k) 的累加
end

az = zeros(1,azord+1); 				% 分母系数向量初始化为零向量
for k = 0:azord					% 依次求各多项式乘积结果
    pln = [1];
    for l = 0:k-1  pln = conv(pln,Nz);   end	% 求Nz的k 次乘积
    pld = [1];
    for l = 0:azord-k-1  pld = conv(pld,Dz);  end		% 求Dz的bzord-k 次乘积
    az = az+aZ(k+1)*conv(pln,pld); 		% 求bZ*Nz^k*Dz^(bzord-k) 的累加
end
az1 = az(1); az = az/az1; bz = bz/az1;		% 将分母多项式首项归一化