MATLAB Program for Response of a 2nd order system

Program:
% transfar function & response of a 2nd order system
% zeta is tha damping factor
% wn is natural frequency
zeta=input('enter zeta ');
wn=input('enter wn in rad/sec ');
num=[wn^2];
den=[1 2*zeta*wn wn^2];
G=tf(num,den)                  %% get transfar function
subplot(2,2,1),step(G)      %% step response response of system
subplot(2,2,2),impulse(G)  %% impulse response of system
subplot(2,2,3),bode(G)     %% bode plot of system
subplot(2,2,4),nyquist(G)  %% nyquist plot of system

Output:

MATLAB Program for Stability of a System Using Routh-Hurwitz Criterion

Stability of any system is an important issue. we have several methods to find out the stability of any system Routh-Hurwitz Criterion is one of them, we can check the stability of system using Routh Matrix.
Routh- Hurwitz Criterion state that 

"The system is stable if and only if all the elements in the first column have same algebaric sign. If all elements are not of the same sign then the number of sign change of the elements  in first column equals the number of rots of characteristic equation in the right half of S-plane"


Using MATLAB can check the stability of any system on the bases of Routh-Hurwitz Criterion with the help of following program:


%% routh hurwitz criteria
clear
clc
%% firstly it is required to get first two row of routh matrix
e=input('enter the coefficients of characteristic equation: ');
disp('-------------------------------------------------------------------------')
l=length(e);
m=mod(l,2);
if m==0
    a=rand(1,(l/2));
    b=rand(1,(l/2));
    for i=1:(l/2)
        a(i)=e((2*i)-1);
        b(i)=e(2*i);
    end
else
    e1=[e 0];
    a=rand(1,((l+1)/2));
    b=[rand(1,((l-1)/2)),0];
    for i=1:((l+1)/2)
        a(i)=e1((2*i)-1);
        b(i)=e1(2*i);
    end
end
%% now we genrate the remaining rows of routh matrix
l1=length(a);
c=zeros(l,l1);
c(1,:)=a;
c(2,:)=b;
for m=3:l
    for n=1:l1-1
        c(m,n)=-(1/c(m-1,1))*det([c((m-2),1) c((m-2),(n+1));c((m-1),1) c((m-1),(n+1))]);
    end
end
disp('the routh matrix:')
disp(c)
%% now we check the stablity of system
if c(:,1)>0
    disp('System is Stable')
else
    disp('System is Unstable');
end

MATLAB Program for Frequency Response Analysis for a Transfer Function

to analysis the frequency response of any transfer function we use FREQS (for analog system or for S-domain) & FREQZ (for digital system or for Z-domain).

Program (for S-domain):

clear
clc
num=input('enter the coefficient of numerator: ');
den=input('enter the coefficient of denominator: ');
freqs(num,den)

Output:

enter the coefficient of numerator: 1
enter the coefficient of denominator: [1 2 8]

Program (for Z-domain): 

clear

clc

num=input('enter the coefficient of numerator: ');

den=input('enter the coefficient of denominator: ');

freqz(num,den)

Output:

enter the coefficient of numerator: [1 1 2]

enter the coefficient of denominator: [2 1 3]

MATLAB Program for Dual pass band Microwave Filter for 6/4 GHz Band

• In satellite communication 6/4 GHz band used in many countries. Commercial communication satellites use a frequency band of 500 MHz bandwidth near 6 GHz for uplink transmission and another 500 MHz bandwidth near near 4 GHz for down link transmission. While an uplink of 5.725 to 7.075 GHz and a downlink of 3.4 to 4.8 GHz is used. For both uplink and downlink band a single dual bandpass filter can be used instead of two band-pass filters.
• To design microwave filter well known classical methods like Butterworth, Chebyshev are used but these methods are applicable on single band filters. To design the multipassband filter some limited analytical techniques have been proposed. These all techniques are based on optimization method. Apart from these techniques another technique to design of microwave filter with arbitrary frequency response also proposed which is based on iterative quasi-Newton algorithm. 
• The methodology proposed in this paper  follow this technique and rests on the translation of analog specification of both bands of 6/4 GHz band into digital specification since in digital domain the design is little bit easy. A designer can also take advantage of these sophisticated and continuously developing digital techniques in digital domain.

Following program can be used to design a dualpassband filter for satellite communication:


%% clear the command window & workspace
clear all
clc
%% define the besic perameter related to filter
% order of numerator % denominator
n=10;
d=10;
% frequency range of passbands & stopbands
f1=0:1/1000:.9406;
f2=0.9406:.0001:.9578;
f3=0.9579:.0001:.9645;
f4=.9646:.0001:.9714;
f5=[.9715:1/1000:1 1];
f=[f1 f2 f3 f4 f5];
edges=[0 f(1278:1280)];
% magnitiude to define passbands & stopbands
m1=[repmat(-40,1,900) -40:1:0];
m2=repmat(0,1,173);
m3=[0:-1:-32 -33:1:0];
m4=repmat(0,1,69);
m5=[0:-1:-27 -28.5 -30];
m=[m1 m2 m3 m4 m5];
%% desigining of filter
% conversion of magnitude from dB to normal value
a=db2mag(m);
% define the weight of error
w=repmat(50,1,1280);
% numerator & denominator of digital filter
[num,den]=iirlpnorm(n,d,f,edges,a,w);
% frequency responce of digital filter
figure(1),freqz(num,den)
% transfer function of digital filter
s=tf(num,den,10^(-9));
% pole zero plot of designed digital filter
figure(2),pzplot(s)
% conversion of digital filter from digital to analog
S=d2c(s,'tustin');
% bode plot of analog filter
figure(3),bode(S)
% numerator & denominator of analog filter
[NUM,DEN]=tfdata(S,'v');
% conversion of negative coefficent of numerator & denominator to make
% filter stable
if NUM<0
    NUM=-1*NUM;
end
if DEN<0
    DEN=-1*NUM;
end
% transfer function of analog filter with new coefficents
S1=tf(NUM,DEN);
% bode plot of analog filter with new coefficents
figure(4),bode(S1)

%% end of designing program

%%output

MATLAB Program for Continuous Time Square Wave

square() function is used to generate the square wave.


%Program to generate a continuous time square wave


close all;
clear all;


a=input('Enter the amplitude of the square wave A = ');
f= input('Enter the frequency of the square wave F = ');
dc=input('Enter the duty cycle of the wave DC = ');
f=f*2*pi;
t=0:.001:1;
y=a*square(f*t,dc);
plot(t,y);
axis([0 1 -2.2 2.2]);


%output


%Enter the amplitude of the square wave A = 2
%Enter the frequency of the square wave F = 10
%Enter the duty cycle of the wave DC = 50

MATLAB Program for Exponentially Damped Sinusoidal Wave

%program to generate an exponentially damped sinusoidal wave


close all;
clear all;


A=input('Enter the amplitude of the sinusoidal wave A = ');
f= input('Enter the frequency of the sinusoidal wave F = ');
phi=input('Enter the phase angle of sinusoid phi = ');
a=input('Enter the attenuation factor a = ');


f=f*2*pi;
t=0:.001:1;
y=A*sin(f*t + phi).*exp(-a*t);


plot(t,y);
axis([0 1 -2.2 2.2]);


%output


%Enter the amplitude of the sinusoidal wave A = 2
%Enter the frequency of the sinusoidal wave F = 10
%Enter the phase angle of sinusoid phi = 0
%Enter the attenuation factor a = 6

MATLAB Program for Linear Convolution of Two sequences

 conv() function is used for finding the linear convolution of sequences. The program is shown below.


%Program to find the linear convolution of two sequences

x1=input('Enter the first sequence x1(n) = ');
t1=input('Enter the starting time of first sequence t1 = ');
x2=input('Enter the second sequence x2(n) = ');
t2=input('Enter the starting time of second sequence t2 = ');
l1=length(x1);
l2=length(x2);
ln=l1+l2-1;

yn=conv(x1,x2);

a=t1+l1-1;
t=t1:a;
subplot(311);
stem(t,x1);
grid on;
xlabel('time--->');
ylabel('amplitude--->');
TITLE('First sequence');

a=t2+l2-1;
t=t2:a;
subplot(312);
stem(t,x2);
grid on;
xlabel('time--->');
ylabel('amplitude--->');
TITLE('Second sequence');

tn=t1+t2;
a=tn+ln-1;
t=tn:a;
subplot(313);
stem(t,yn);
grid on;
xlabel('time--->');
ylabel('amplitude--->');
TITLE('Convolved output');

%output

%Enter the first sequence x1(n) = [1 2 6 2 3 1 4]
%Enter the starting time of first sequence t1 = -3
%Enter the second sequence x2(n) = [3 1 4 5 2]
%Enter the starting time of second sequence t2 = -1

MATLAB Program for FIR Bandpass Filter

function fir1() is used for the design of FIR filter.You can refer the MATLAB product help for the syntax of these functions. The program for an FIR bandpass filter is given below.

% Program to design FIR bandpass filter .

close all;
clear all;

fp=input('Enter the start and stop frequency of the pass band');
f=input(' Enter the sampling frequency');
n=input(' Enter the order of the filter');

% Normalizing the frequencies

wp=(2/f).*fp;

%Calculation of filter coefficients

b=fir1(n,wp);

%Plotting the filter response

freqz(b,1,500,f);
TITLE('Magnitude and Phase response');


%output
%Enter the start and stop frequency of the pass band[1000 2000]
%Enter the sampling frequency 5000
%Enter the order of the filter 50

MATLAB Program for IIR Chebyschev Type - II Filter

The function cheb2ord() is used to find the filter order and center frequency and cheby2() is used to find the filter coefficients.A low pass filter is implimented.


% Program to design IIR Chebyschev type - II filter


clear all;
close all;


fp=input('Enter the pass band frequency fp   = ');
fs=input('Enter the stop band frequency fs   = ');
rp=input('Enter the pass band attenuation rp = ');
rs=input('Enter the stop band attenuation rs = ');
f=input ('Enter the sampling frequency f     = ');


wp=2*fp/f;
ws=2*fs/f;


[n,wn]=cheb2ord(wp,ws,rp,rs);


[b,a]=cheby2(n,rs,wn,'low');


freqz(b,a,500,f);
TITLE ('Magnitude and phase respose of the IIR Chebyschev type - II filter');


%output
%Enter the pass band frequency fp   = 1000
%Enter the stop band frequency fs   = 1200
%Enter the pass band attenuation rp = .2
%Enter the stop band attenuation rs = 45
%Enter the sampling frequency f     = 3000

MATLAB Program for IIR Chebyschev Type - I Filter

The function cheb1ord() is used to find the filter order and center frequency and cheby1() is used to find the filter coefficients.A low pass filter is implimented.

% Program to design IIR Chebyschev type - I filter


clear all;
close all;


fp=input('Enter the pass band frequency fp   = ');
fs=input('Enter the stop band frequency fs   = ');
rp=input('Enter the pass band attenuation rp = ');
rs=input('Enter the stop band attenuation rs = ');
f=input ('Enter the sampling frequency f     = ');


wp=2*fp/f;
ws=2*fs/f;


[n,wn]=cheb1ord(wp,ws,rp,rs);


[b,a]=cheby1(n,rp,wn,'low');


freqz(b,a,500,f);
TITLE ('Magnitude and phase respose of the IIR Chebyschev type - I filter');


%output
%Enter the pass band frequency fp   = 1000
%Enter the stop band frequency fs   = 1200
%Enter the pass band attenuation rp = .2
%Enter the stop band attenuation rs = 45
%Enter the sampling frequency f     = 3000

MATLAB Program for IIR Butterworth Filter

The function buttord(), is used to find the filter order and center frequency . Function butter() is used to find the filter coefficients.A low pass filter is implemented.


%Program to design an IIR Butterworth filter


clear all;
close all;
fp=input('Enter the pass band frequency fp   = ');
fs=input('Enter the stop band frequency fs   = ');
rp=input('Enter the pass band attenuation rp = ');
rs=input('Enter the stop band attenuation rs = ');
f=input ('Enter the sampling frequency f     = ');


wp=2*fp/f;
ws=2*fs/f;


[n,wn]=buttord(wp,ws,rp,rs);


[b,a]=butter(n,wn,'low');


freqz(b,a,500,f);
TITLE ('Magnitude and phase respose of the IIR butterworth filter');


%output
%Enter the pass band frequency fp   = 1000
%Enter the stop band frequency fs   = 1200
%Enter the pass band attenuation rp = .2
%Enter the stop band attenuation rs = 45
%Enter the sampling frequency f     = 3000

MATLAB Program for FIR Filter Using Window

Popular window coefficients.

1. hann() - for hanning window
2. hamming() - for hamming window
3. blackman() - for blackman window
4. kaiser() - for kaiser window

% Program to design a FIR filter using windows.

close all;

clear all;

fp=input('Enter the pass band frequency');

fs=input('Enter the stop band frequency');

rp=input(' Enter the pass band attenuation');

rs=input('Enter the stop band attenuation');

f=input(' Enter the sampling frequency');

%Calculating filter order

num=-20*log10(sqrt(rp*rs))-13;

dem=14.6*(fs-fp)/f;

n=ceil(num/dem);

n=abs(n);

% Normalizing the frequencies

wp=2*fp/f;

ws=2*fs/f;

wn=(ws+wp)/2;

%Adjusting the filter order. The order of window must be an odd number 

%and the order of filter must be one less than that of the window 

if (rem(n,2)==0)

    m=n+1;

else

        m=n;

        n=n-1;

end

%Window sequence calculation

w=hann(m);

%Calculation of filter coefficients

b=fir1(n,wn,'low',w);

%Plotting the filter response

freqz(b,1,500,3000);

TITLE('Magnitude and Phase response');

Output:

%output

%Enter the pass band frequency1000

%Enter the stop band frequency1200

%Enter the pass band attenuation.2

%Enter the stop band attenuation45

%Enter the sampling frequency3000

You can change this lowpass filter to high pass filter by changing the option 'low' to 

'high' in the fir1() function. The output is shown below.

%output

%Enter the pass band frequency1200

%Enter the stop band frequency1000

%Enter the pass band attenuation.2

%Enter the stop band attenuation45

%Enter the sampling frequency3000