function index=autodesign(num,den,wm,kc,op)
%==========================================================================
%自動(dòng)控制原理課程設(shè)計(jì)實(shí)驗(yàn)程序。
%輸入變量num,den分別為開(kāi)環(huán)環(huán)傳遞函數(shù)展開(kāi)式中分子與分母系數(shù)
%校正采用頻率法超前校正，wm為期望截止頻率，當(dāng)不實(shí)行校正時(shí)wm置0
%kc為根軌跡增益，kc取不同的值可得到不同增益校正后的傳遞函數(shù)，若不用則置1
%op為操作項(xiàng)，op為0時(shí)，畫(huà)閉環(huán)階躍響應(yīng)曲線，op為1時(shí)，畫(huà)波德圖，op為2時(shí)畫(huà)根軌跡，op為3時(shí)對(duì)原系統(tǒng)進(jìn)行校正，且畫(huà)校正后階躍響應(yīng)曲線，此時(shí)wm不能為0
%輸出為時(shí)域和頻域能指標(biāo)，格式為index=[td tr tp ts siagama r h wc wx]
%當(dāng)輸出index矩陣有兩行時(shí)，第一行為校正前的指標(biāo)，等二行為校正后的指標(biāo)
%當(dāng)op不為3時(shí)，則不進(jìn)行校正操作，index只有一行，為原系統(tǒng)性能指標(biāo)
%==========================================================================
%==========================================================================
if op~=3
    wm=0;
end
go=tf(num,den);%求開(kāi)環(huán)傳遞函數(shù)
go=kc*go;
[index,mag,w]=ide(go); %求當(dāng)前系統(tǒng)各項(xiàng)時(shí)域和頻域指標(biāo)
%畫(huà)階躍響應(yīng)曲線


%頻率法校正
if wm~=0
    [wmmin,loc]=min(abs(w-wm));
    wm=w(loc);
    a=(1/mag(loc))^2;
    T=1/(wm*sqrt(a));
    gc=tf([a*T 1],[T 1]);
    gcgo=gc*go;
    index=[index;ide(gcgo)];%求校正后系統(tǒng)的各項(xiàng)指標(biāo)
end

switch op
    case 0        
        go=go/(1+go);
        [y,t]=step(go,[0:0.1:40]);  %從0s,到100s,以0.1s為精度做階躍響應(yīng)實(shí)驗(yàn)
        plot(y);hold on;plot([0:400],ones(1,401),':r'); %畫(huà)階躍響應(yīng)曲線
    case 1
        bode(go);grid on;%畫(huà)波得圖
    case 2
        rlocus(go); %做根軌跡
    case 3
        gcgo=gcgo/(1+gcgo);
        [y,t]=step(gcgo,[0:0.1:40]);  %從0s,到100s,以0.1s為精度做階躍響應(yīng)實(shí)驗(yàn)
        plot(y);hold on;plot([0:400],ones(1,401),':r'); %畫(huà)階躍響應(yīng)曲線
end
        
function [index,mag,w]=ide(sys)%求該傳遞函數(shù)性能指標(biāo)（頻域及時(shí)域）
sysp=sys/(1+sys);
[y,t]=step(sysp,[0:0.1:100]);  %從0s,到100s,以0.1s為精度做階躍響應(yīng)實(shí)驗(yàn)             
len=length(y);
yend=y(len);                %取終值
ymax=max(y);                %取最大值（峰值）
sigama=(ymax-yend)/yend;    %計(jì)算超調(diào)量

%計(jì)算調(diào)節(jié)時(shí)間---------------------------------------
for ii=len:-1:1
    if abs(y(ii)-yend)>0.05*yend   
        ts=t(ii);               
        break
    end
end

for ii=2:len                      %計(jì)算峰值時(shí)間
    if y(ii)-y(ii-1)<=0
        tp=t(ii);
        break
    end
end

%計(jì)算上升時(shí)間-----------------------------------------
for ii=2:len           
    if y(ii)>=0.1*yend & y(ii-1)<0.1*yend
        tr1=t(ii);
    end
    if y(ii)>0.9*yend & y(ii-1)<=0.9*yend
        tr2=t(ii-1);
        break
    end
end
tr=tr2-tr1;                              
%計(jì)算延遲時(shí)間----------------------------------------
for ii=2:len
    if y(ii)>=0.5*yend & y(ii-1)<0.5*yend  
        td=t(ii); 
        break
    end
end
%頻域分析--------------------------------------------         
[mag phase w]=bode(sys);
[magmin,loc]=min(abs(mag-1));
wc=w(loc);            %計(jì)算截止頻率
r=180+phase(loc);     %計(jì)算相角裕度
[phasemin,loc]=min(abs(phase+180));
wx=w(loc);           %計(jì)算穿越頻率
h=1/abs(mag(loc));  %計(jì)算幅值裕度
index=[td tr tp ts sigama r h wc wx];