阅读量:0
相关代码打开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 end function [h, w] = freqziir(b, a, count, whole) if whole lastpoint = 2 * pi; else lastpoint = pi; end w = [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; end function [h, w] = freqzfir(b) a = [1]; count = 512; n = length(b); whole = false; if whole lastpoint = 2 * pi; else lastpoint = pi; end w = [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