evaltfp.m

高阶sigma-delta调制器设计matlab工具包, 半波带滤波器设计工具包

function H=evalTFP(Hs,Hz,f)
% H=evalTFP(Hs,Hz,f)
% Compute the value of a transfer function product Hs*Hz at a frequency f,
% where Hs is a cts-time TF and Hz is a discrete-time TF.
% Both Hs and Hz are SISO zpk objects.
% This function attempts to cancel poles in Hs with zeros in Hz.

szeros = Hs.z{1};
spoles = Hs.p{1};
zzeros = Hz.z{1};
zpoles = Hz.p{1};

H = zeros(size(f));
w = 2*pi*f;	s = j*w;	z=exp(s);
for i=1:length(f)
    wi = w(i);	si = s(i);	zi = z(i);
    if isempty(spoles)
        cancel = 0;
    else
        cancel = abs(si-spoles)<1e-5;
    end
    if ~cancel
    	% wi is far from a pole, so just use the product Hs*Hz
    	H(i) = evalTF(Hs,si) * evalTF(Hz,zi);
    else
    	% cancel pole(s) of Hs with corresponding zero(s) of Hz
		cancelz = abs(zi-zzeros)<1e-5;
		if sum(cancelz) > sum(cancel)
		    H(i) = 0;
		elseif sum(cancelz) < sum(cancel)
	    	    H(i) = Inf;
		else
		    H(i) =  evalRPoly(szeros,si,Hs.k) * ...
			    zi^sum(cancel) * evalRPoly(zzeros(~cancelz),zi,Hz.k) / ... 
			    (evalRPoly(spoles(~cancel),si,1)*evalRPoly(zpoles,zi,1));
		end
    end
end