why is my phase spectrum calculated by hand not agreeing with matlab's freqz?

I am trying to replicate the Matlab function freqz. To calculate the phase of the frequency response of my filter I am implementing the equation for phase:

$$\phi(e^{j\omega}) = \tan{}^{-1}\left\{\frac{\Im{}[H(z)]}{\Re{[H(z)]}}\right\}_{z=e^{j\omega}}$$

and then unwrapping it by hand.

Here is the code I have:

%%
[b,a] = cheby2(4, 40, .5, 'low');
[Hz fVec] = freqz(b,a);
resolution = length(fVec);
% normalized frequency scale to plot against
scale = linspace(0, 1-1/resolution, resolution);
% make a vector of complex frequencies from 0 to (almost) pi
omega = exp(i*scale*pi);
% find the poles and zeros of the filter
zeros = roots(b);
poles = roots(a);
% compute magnitude and phase at each frequency
for n=1:resolution
num = 1;
for m=1:length(zeros)
% multiply the polynomial (1-zq0)...(1-zqn)
num = num*(1-omega(n)*zeros(m));
end
den = 1;
for m=1:length(poles)
% multiply the polynomial (1-zp0)...(1-zpn)
den = den*(1-omega(n)*poles(m));
end
phase(n) = atan(imag(num/den)/real(num/den));
mags(n) = 20*log10(abs(num/den));
end
% unwrap by hand
for n=2:length(phase)
while abs(phase(n-1)-phase(n)) > (pi/2)
if phase(n-1)>phase(n)
phase(n) = phase(n) + pi;
else
phase(n) = phase(n) - pi;
end
end
end
% plot and compare
figure;
subplot(2,1,1)
plot(scale, mags);
subplot(2,1,2)
plot(scale, phase);
figure;
freqz(b,a);


Now I expect my plot to be more or less identical to that plotted by freqz (some plot annotations and phase expressed in radians rather than degrees aside). But it is not:

The magnitude spectrum is indeed identical apart from a gain factor, but the phase quite different. Mine doesn't have those kinks in it, but that is because of my unwrapping. However, the sign is reversal is a mystery to me. Can anyone explain the reason for this discrepancy?

1 Answer

A few things to note here

1. The polynomials are written in $$z^{-1}$$ and not in $$z$$. You need to take the inverse, it should be omega = 1./exp(1i*scale*pi); or, simpler omega = exp(-1i*scale*pi); . Frankly, 'omega' is not a good variable name for this. It should be "z" or "zToTheMinus1" instead.
2. You loose the gain in the process so your overall scaling is way off
3. atan() can only resolve angles from $$-\pi$$ to $$+\pi$$, use atan2() instead
4. You use degrees on one graph and radians on the other.
• 1 is what I had missed. Thank you. In these days of not having colleagues in the room to ask to look at things like this it is very great to have resources like this. I think I mentioned points 2 and 4 in the question. I am using atan and not atan2 because that is equivalent to the function that I am using in non-Matlab land in this application. I believe the wrapping is correct, however. Aug 18, 2021 at 6:19
• Cool. atan2() is easy enough to write if you have access to atan() and it's much safer to use. Your unwrapping only works if your first phase calculation is correct, so you lucked out here. atan2() also handles $X_{real}= 0$ correctly, whereas trying to use atan() will create a "divide by zero" error Aug 18, 2021 at 11:22