0
$\begingroup$

I am making some homemade tools on Matlab. I made function that plot the frequency response of a discrete transfer function (just like freqz does). When I compare my result with the result from freqz, I get something very similar, but not exactly the same.

Can somebody see why ? I am a little confused here.

Observations

  • I first thought it was related to frequency warping, but when I applied warping on my frequencies, I got a similar result.
  • We also see that the phase does not wrap correctly, it passes from -179.3 to 180.7 (instead of a perfect -180 - 180);

Graph enter image description here

Code

function [h,f] = py_freqz(b,a,fs,n)
%make sur both polynomials are same size
if (length(a) > length(b))
    b = [b, zeros(1,max(0,length(a)-1))]; 
elseif (length(b) > length(a))
    a = [a, zeros(1,max(0,length(b)-1))]; 
end
p = roots(a);
z = roots(b);

precision = fs/n;
w = linspace(0,pi-precision/2,n);
f = w/pi*fs/2;
h = ones(1,n);
phases = zeros(1,n);
unitcircle = exp(1i*w);
%Zeros
for k=1:length(z)
    v = unitcircle-z(k);
    h = h .* abs(v);
    phases = phases + angle(v);
end
%Poles
for k=1:length(p)
    v = unitcircle-p(k);
    h = h ./ abs(v);
    phases = phases - angle(v);
end

h = h * b(1);   % Input Gain

subplot(2,1,1);
xlabel('Frequency [Hz]');
ylabel('Magnitude');
plot(f,10*log(h));
grid on

subplot(2,1,2);
xlabel('Frequency [Hz]');
ylabel('Phase');
plot(f,phases/pi*180);
grid on
end

I generated the above graph using

py_freqz([0.09247185865855016 2*0.09247185865855016 0.09247185865855016 ], [1 -1.5668683171334845 0.9367557517676852], 1000,1000);
$\endgroup$
3
$\begingroup$

Several corrections:

  1. This code does not make sense:
precision = fs/n;
w = linspace(0,pi-precision/2,n);
f = w/pi*fs/2;

Your precision should be the resolution in frequency domain (I think). But, now consider n=fs. Then precision=1 (which is fine), but w ranges only from (0...pi-0.5), i.e. you do not span the entire unit circle. The corrected code is (you actually dont need the precision variable):

w = linspace(0, pi, n+1); w = w(1:end-1);
f = w/pi*fs/2;

The trick is to create n+1 points in the interval [0,pi] and the remove the last one. This is much simpler than calculating the offset on your own.

  1. You cannot add up the phases from each frequency zero and pole directly. Instead, you would need to do a modulo-addition (i.e. wrap around the phase). However, I suggest to simply track the phase+magnitude in a single variable h_all:
    unitcircle = exp(1i*w);
    h_all = b(1)*ones(1,n);
    %Zeros
    for k=1:length(z)
        v = unitcircle-z(k);
        h_all = h_all .* v;
    end
    %Poles
    for k=1:length(p)
        v = unitcircle-p(k);
        h_all = h_all ./ v;
    end

And finally you use this to extract phase and magnitude.

h = abs(h_all);
phases = angle(h_all);

Here's the full code for the function:

function [h_all,f] = py_freqz(b,a,fs,n)
%make sur both polynomials are same size
if (length(a) > length(b))
    b = [b, zeros(1,max(0,length(a)-1))]; 
elseif (length(b) > length(a))
    a = [a, zeros(1,max(0,length(b)-1))]; 
end
p = roots(a);
z = roots(b);

w = linspace(0, pi, n+1); w = w(1:end-1);
f = w/pi*fs/2;
unitcircle = exp(1i*w);
h_all = b(1)*ones(1,n);
%Zeros
for k=1:length(z)
    v = unitcircle-z(k);
    h_all = h_all .* v;
end
%Poles
for k=1:length(p)
    v = unitcircle-p(k);
    h_all = h_all ./ v;
end

subplot(2,1,1);
xlabel('Frequency [Hz]');
ylabel('Magnitude');
plot(f,10*log(abs(h_all)));
grid on

subplot(2,1,2);
xlabel('Frequency [Hz]');
ylabel('Phase');
plot(f,angle(h_all));
grid on
end

And the test script

b = [0.09247185865855016 2*0.09247185865855016 0.09247185865855016 ];
a = [1 -1.5668683171334845 0.9367557517676852];

[h, w] = freqz(b, a);
[h2_all, w2] = py_freqz(b, a, 1,8192);

%close all;
figure(1)
subplot(211);
plot(w/(2*pi), 20*log10(abs(h)), 'b-x'); 
hold on;
plot(w2, 20*log10(abs(h2_all)), 'r-');
hold off;

subplot(212);
hold off;
plot(w/(2*pi), angle(h), 'b-x');
hold on;
plot(w2, angle(h2_all), 'r-');

And the program output:

enter image description here

$\endgroup$
  • $\begingroup$ Thank you for this. As for the first piece of code, I think it used to make sense, but I am pretty sure I messed up the code right before posting in on Stack Exchange. Working with complex number and ploting abs() and angle() seems to have fixed it. $\endgroup$ – Pier-Yves Lessard Jan 19 '17 at 2:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.