# 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. 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.
3. atan() can only resolve angles from $$-\pi$$ to $$+\pi$$, use atan2() instead
• 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