num=1000;den=conv([1,0],conv([0.1,1],[0.001,1]));
G=tf(num,den);                %未校正系統(tǒng)的開環(huán)傳遞函數(shù)
[Gm,Pm,Wcg,Wcp]=margin(G);  %1) 未校正系統(tǒng)的頻域響應(yīng)參數(shù)，計(jì)算需要相角欲度Pm(γ)
w=0.1:0.1:10000;              %確定頻率的取值范圍和頻率采樣的間隔值
[mag,phase]=bode(G,w);magdb=20*log10(mag);  %對數(shù)幅頻、相頻特性向量，計(jì)算需要magdb
%(對數(shù)幅頻響應(yīng)值)
phim1=45;deta=8;              % 2）(4.1.9)式設(shè)置系統(tǒng)設(shè)計(jì)參數(shù)phim1( )；deta( )
phim=phim1-Pm+deta;                        % (4.1.9)式相位超前角φ 
bita=(1-sin(phim*pi/180))/(1+sin(phim*pi/180));        %求出β值
n=find(magdb+10*log10(1/bita)<=0.0001);       %找出滿足magdb+10*log10(1/bita)<=0.0001
%式的magdb向量所有的下標(biāo)值
wc=n(1);      %通常magdb(1)+10*log10(1/bita)>0.0001;取 的第1項(xiàng)為 (比實(shí)際頻率值大
%10倍)，這是因?yàn)閣=0.1:0.1:10000，而下標(biāo)向量 的值=1:1:100000
w1=(wc/10)*sqrt(bita);w2=(wc/10)/sqrt(bita);               % 3)  (4.1.11)和(4.1.12)式
numc=[1/w1,1];denc=[1/w2,1];                           % (4.1.13)式，取K=1 
Gc=tf(numc,denc);                                     %校正器的傳遞函數(shù)                       
GmdB=20*log10(Gm);                         %4) 下面程序校驗(yàn)系統(tǒng)校正后的系統(tǒng)指標(biāo)
GcG=Gc*G;[Gmc,Pmc,wcgc,wcpc]=margin(GcG);        %GcG是校正后系統(tǒng)的開環(huán)傳遞函數(shù)
GmcdB=20*log10(Gmc);                       
disp('未校正系統(tǒng)的開環(huán)傳遞函數(shù)和頻域響應(yīng)參數(shù)：h，γ，wc')
G,[GmdB,Pm,Wcp],
disp('校正裝置傳遞函數(shù)和校正后系統(tǒng)開環(huán)傳遞函數(shù)')
Gc,GcG,
disp('校正后系統(tǒng)的頻域響應(yīng)參數(shù)：h，γ，wc')
[GmcdB,Pmc,wcpc],
disp('校正裝置的參數(shù)T和β值：T，β')
T=1/w1;[T,bita],
bode(G,GcG);figure(2);margin(GcG)