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.
- 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);
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);