function faza;
b=[1 3 3 1];
a=[1 -3 4 -2 1];
[h,w]=freqs(b,a);
phi=angle(h);
phi=unwrap(phi);
% plot(w,phi*180/pi);
% plot(w,abs(h));
freqs(b,a);
grid on