% frequency response of LTI systems
%
% @Kefu Xue, Ph.D., Apr. 2001
%
echo on;
% frequency response and Bode plot
%
% H(s)=20s(s+100)/((s+2)(s+10))
%
wr=0.1;% frequency resolution
w=0.01:wr:10000; %0.01 to 10000 radian/second
wlog=log10(w); %log base 10 frequency axis
s=j*w;
H=20*s.*(s+100)./((s+2).*(s+10));
M=abs(H);
Phi=angle(H);
subplot(2,1,1);plot(w,M);title('magnitude frequency response');
xlabel('radian/second');
subplot(2,1,2);plot(w,Phi);title('phase frquency response');
ylabel('radian');xlabel('radian/second');
pause
%
% using dB and log10 scale in frequency
%
Mdb=20*log10(M);% in dB
Adeg=180*Phi/pi;% in degree
figure;subplot(2,1,1);plot(wlog,Mdb);title('magnitude frequency response');
ylabel('dB');xlabel('10^( ^) radian/second');
subplot(2,1,2);plot(wlog,Adeg);title('phase frquency response');
ylabel('degree');xlabel('10^( ^) radian/second');
pause
%
% H(s)=s^2/((s+1)(s^2+4s+16))
%
echo off;
H=s.^2./((s+1).*(s.^2+4*s+16));
M=abs(H);
Phi=angle(H);
Mdb=20*log10(M);% in dB
Adeg=180*Phi/pi;% in degree
figure;subplot(2,1,1);plot(wlog,Mdb);title('magnitude frequency response');
ylabel('dB');xlabel('10^( ^) radian/second');
subplot(2,1,2);plot(wlog,Adeg);title('phase frquency response');
ylabel('degree');xlabel('10^( ^) radian/second');
pause
echo on;
%
%Sinusoidal Steady State Response 
%
% H(s)=20/(s+20)is system transfer function
% input signal x(t)=cos(4pi*t)+sin(50pi*t)
%
w0=4*pi; 
T0=2*pi/w0; %find period of signal x(t)
deltaT=0.01;
t=0:deltaT:2*T0; % for 2 periods of x(t)
xt=cos(4*pi*t)+sin(50*pi*t); %input signal
%
egvalu1=20/(j*4*pi+20); %eigenvalue for w=4pi
egvalu2=20/(j*50*pi+20); %eigenvalue for w=50pi
%
ysst=abs(egvalu1)*cos(4*pi*t+angle(egvalu1))+abs(egvalu2)*sin(50*pi*t+angle(egvalu2));
%
figure;subplot(2,1,1);plot(t,xt,'r--',t,ysst,'b');legend('x(t)','yss(t)');
f=1:0.1:30; %frequency 1 to 30 Hz
s=j*2*pi*f;
H=20./(s+20);
subplot(2,1,2);plot(f,abs(H));title('magnitude frequency response');
pause
%
% H(s)=(-s+10)/(s+10)
% input signal x(t)=cos(4pi*t)+cos(12pi*t)+sin(32pi*t)
echo off
xt=cos(4*pi*t)+cos(12*pi*t)+sin(32*pi*t); %input signal
%
egvalu1=(-j*4*pi+10)/(j*4*pi+10); %eigenvalue for w=4pi
egvalu2=(-j*12*pi+10)/(j*12*pi+10); %eigenvalue for w=12pi
egvalu3=(-j*32*pi+10)/(j*32*pi+10); %eigenvalue for w=32pi
%
% steady state responses of signal components
%
x1=cos(4*pi*t);
yss1=abs(egvalu1)*cos(4*pi*t+angle(egvalu1));
x2=cos(12*pi*t);
yss2=abs(egvalu2)*cos(12*pi*t+angle(egvalu2));
x3=sin(32*pi*t);
yss3=abs(egvalu3)*sin(32*pi*t+angle(egvalu3));
figure;subplot(2,2,1);plot(t,x1,'r--',t,yss1,'b');
subplot(2,2,2);plot(t,x2,'r--',t,yss2,'b');
subplot(2,2,3);plot(t,x3,'r--',t,yss3,'b');
ysst=yss1+yss2+yss3;
%
subplot(2,2,4);plot(t,xt,'r--',t,ysst,'b');legend('x(t)','yss(t)');
pause
H=(-s+10)./(s+10);
figure;subplot(2,1,1);plot(f,abs(H));title('magnitude frequency response');
subplot(2,1,2);plot(f,angle(H));title('phase frequency response');
%
pause
close all
