An Assignment on Bode Plots

This post is something to do with an assignment I made for someone else. It is related to subject Modern Control Systems that I studied in 1999 for my BE Avionics. The instructor used to give assignments for the course. For one assignment I got the remarks “No Copying Next Time” by the instructor. I was stunned. The truth was the other way round. I got the wind of the news that someone from the class had copied my assignment from my account in Windows NT Workstation with the assistance of the computer lab staff. USB flash drives had not been invented yet. For next assignment I made about 75% of the assignment and completed the rest on final day. I was on the top of the class in the subject. The last assignment was bit challenging and final exams were close. I decided not to make the assignment. Whole class did not submit. There was our course mate who was struggling in the course. The instructor told him that if you submit last assignment you can be given C. He came to me and requested that I should make assignment for him. I agreed. I made this assignment for him. He got C for the course. As far as i can remember that there were more than 10 Ds in the course. He told me that his Pak No in fruit shop is open for me. I told him that I do not need that and asked him to give me just one hard copy of the assignment. I am sharing it as i considered it to be my best piece of MATLAB code I wrote during my BE. I only have the function definition file. With the help of this function some example graph were drawn. The main theme of the assignment was to develop algorithm for drawing approximate bode plot and then compare it with actual bode plot.

function [errg,errp]=bodeplots(nbr,dbr,g,w)
c=g;
lnb=length(nbr);
ldb=length(dbr);
if dbr~='k'
for n=1:ldb;
if dbr(n)~=0
c=c*dbr(n);
end
end
end
if nbr~='k'
for n=1:lnb;
if (nbr(n)~=0)
c=(1/nbr(n))*c;
end
end
end
if nbr~='k'
num=c*poly(-1*nbr);
else
num=c;
end
if dbr~='k'
den=poly(-1*dbr);
else
den=1;
end
[m1,p1,w]=bode(num,den,w);
db1=20*log10(m1);
db2=20*log10(g)*ones(size(w));
if nbr~='k'
for n=1:lnb;
if nbr(n)==0
temp=20*log10(w);
db2=db2+temp;
else
temp1=20*log10(w/abs(nbr(n)));
temp2=(temp1>=0).*temp1;
db2=db2+temp2;
end
end
end
if dbr~='k'
for n=1:ldb;
if dbr(n)==0
temp=20*log10(w);
db2=db2-temp;
else
temp1=20*log10(w/abs(dbr(n)));
temp2=(temp1>=0).*temp1;
db2=db2-temp2;
end
end
end
dberr=abs(db1-db2);
errg=max(dberr);
figure(1)
subplot(2,1,1),semilogx(w,db1,'-k',w,db2,':k');
legend('Actual','Asympotitic');
title('Bode Magnitude Plot'),ylabel('Gain(db)');
subplot(2,1,2),semilogx(w,dberr,'-k');
xlabel('Freq(rad/s)'),ylabel('Error(db)')
pause
p2=zeros(size(w));
if nbr~='k'
for n=1:lnb;
if nbr(n)==0
temp=90*ones(size(w));
p2=p2+temp;
else
temp=45*log10(w/(0.1*abs(nbr(n))));
temp1=(temp>=0).*temp;
temp2=(temp1<=90).*temp1+(temp1>90)*90;
p2=p2+(nbr(n)/abs(nbr(n)))*temp2;
end
end
end
if dbr~='k'
for n=1:ldb;
if dbr(n)==0
temp=90*ones(size(w));
p2=p2-temp;
else
temp=45*log10(w/(0.1*abs(dbr(n))));
temp1=(temp>=0).*temp;
temp2=(temp1<=90).*temp1+(temp1>90)*90;
p2=p2-(dbr(n)/abs(dbr(n)))*temp2;
end
end
end
perr=abs(p1-p2);
errp=max(perr);
figure(2)
subplot(2,1,1),semilogx(w,p1,'-k',w,p2,':k');
legend('Actual','Asympotitic');
title('Bode Phase Plot'),ylabel('Angle(deg)');
subplot(2,1,2),semilogx(w,perr,'-k');
xlabel('Freq(rad/s)'),ylabel('Error(deg)')

The code is provided as is. It is now difficult for me to explain what i wrote about 12 years ago. A pdf file can be downloaded from here.