相关代码打开matlab源码也可以看到,这里做了简单实现,与源码并不完全一样。
实现代码
[h2 w2] = freqzfir(data);
[h1 w1] = freqz(data);
h2=h2';
h12 = [h1, h2];[h4 w4] = freqziir(b,a, 2001,true);
[h3 w3] = freqz(b,a, w4', 'whole');
h4 = h4';
h34 = h3-h4;function y = mpolyval(x, c) y = zeros(size(x));y(:) = c(1);for i = 2:length(c)y = x .* y + c(i);end
endfunction [h, w] = freqziir(b, a, count, whole)if wholelastpoint = 2 * pi;elselastpoint = pi;endw = [0:count-1]./count.*lastpoint;
% w = w';c = 0 - 1i;zm1 = exp(c * w); h1 = mpolyval(zm1, a);h2 = mpolyval(zm1, b);h = h2./h1;
endfunction [h, w] = freqzfir(b)a = [1];count = 512;n = length(b);whole = false;if wholelastpoint = 2 * pi;elselastpoint = pi;endw = [0:count-1]./count.*lastpoint;
% w = w';c = 0 - 1i;zm1 = exp(c * w); mh2 = mpolyval(zm1, b); mh1 = exp(-1i*w*(n-1));mh = mh2./mh1;h = mh;
end