Matlab freqz 代码简单实现

相关代码打开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