close all; 
clear all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%系統離散化及初始值設置%%%%%%%%%%%%%%%%%%%%%
F0=1;
F1=350;
N=20*F1;
Fs=10*F1;
%sys=tf([-5],[2.0e-5 2e-9 1]);
sys=tf([3.355e7],[1 1.504e3 5.394e5 3.291e7]);
sysd=c2d(sys,1/5000,'tustin');
[num,den]=tfdata(sysd,'v');
[h0,f0]=freqz(num,den,N/2,Fs);
mag=abs(h0);
ph=angle(h0);
ph=ph*180/pi;
ph=ph*180/pi;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%輸入信號%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  A=0.3;%頻率F1信號的幅度
    f0=1
    f1=700;
    N1=20*f1;
    N=30*f1;
    Fs=10*f1;
    dfs=0.5;
    t1=[0:1/Fs:N1/Fs]; %采樣時刻
    t=[0:1/Fs:N/Fs]; 
    F=([1:N]-1)*Fs/N; %換算成實際的頻率值
    F=F(1:N/2);%取N/2個實際頻率點
    df=4*dfs;%頻率間隔
    S=0;
    for i=0:1:((f1-f0)/df)
     f=f0+df*i;
     S=A*cos(2*pi*f*t1)+S;  %組合正弦波
    end;
    x=8*chirp(t1,f0,N/Fs,f1);%chirp
    for i=N1:N+1 %補零消除柵欄效應
    S(i)=0;
    x(i)=0;
    end;

 [Am,Pm,A1m,P1m,A0m,P0m,ym]=qiuzhi(S,N,num,den);
 [Axm,Pxm,A1xm,P1xm,A0xm,P0xm,y1m]=qiuzhi(x,N,num,den);
 k1=round(f0/(Fs/N));
 k2=round(f1/(Fs/N));
pr=1;
for i=k1:k2
   A(pr)=Am(i);
   P(pr)=Pm(i);
   Ax(pr)=Axm(i);
   Px(pr)=Pxm(i);
   ff(pr)= F(i);
   pr=pr+1;
end;