% Name: frequs.m
% ===============================================================
clear,close all
num=input('input the numerator polynomial coefficient vector =:');
den=input('input the denominator polynomial coefficient vector =:');
reg=input('input the frequency region w=');
w=-3*reg:0.01*reg:3*reg;

p=roots(den);                       %求系統極點
q=roots(num);                       %求系統零點
p=p';					          %將極點列向量轉置為行向量
q=q';					          %將零點列向量轉置為行向量
x=max(abs([p q]));				  %確定縱坐標范圍
x=x+1;
y=x;						      %確定橫坐標范圍
subplot(2,2,1);
axis([-x x -y y]);		          %確定坐標軸顯示范圍
axis('square')
plot([-x x],[0 0]);hold on; 	       %畫橫坐標軸
plot([0 0],[-y y]);hold on;		       %畫縱坐標軸
plot(real(p),imag(p),'x');hold on;     %畫極點
plot(real(q),imag(q),'o');hold on;     %畫零點
grid on;
title('zero-pole and frequency response');
ylabel('zero-pole diagram');

[H,w]=freqs(num,den,w);
subplot(2,2,2)
plot(w,abs(H));
grid on;title('magnitude response');
axis([min(w),max(w),min(abs(H)),max(abs(H))])
subplot(2,2,4)
plot(w,angle(H));grid on;
xlabel('frequency radians/s ');title('phase response');
axis([min(w),max(w),min(angle(H)),max(angle(H))])
subplot(223)
t=0:0.01:4/0.1;
h=impulse(num,den,t);
plot(t,h)
title('Impulse response h(t)')
xlabel('Time t')