function [bZ,aZ]=zmapping(bz,az,NZ,DZ)
% 從z域到Z域的頻率變換
% 實現:
% b(Z)/a(Z)=b(z)/a(z),when z=N(Z)/D(Z)

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 i=0:k-1
        pln=conv(pln,NZ);
    end
    pld=[1];
    for i=0:bzord-k-1
        pld=conv(pld,DZ);
    end
    bZ=bZ+bz(k+1)*conv(pln,pld);
end
aZ=zeros(1,azord+1);
for k=0:azord
    pln=[1];
    for i=0:k-1
        pln=conv(pln,NZ);
    end
    pld=[1];
    for i=0:bzord-k-1
        pld=conv(pld,DZ);
    end
    aZ=aZ+az(k+1)*conv(pln,pld);
end    
aZ1=aZ(1);
aZ=aZ/aZ1;
bZ=bZ/aZ1;