%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%   組合正弦波和chirp法測頻率特性MATLAB仿真    %%%%%%%%%
%%%%%%     虛擬儀器頻率特性測試方法仿真            %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; 
clear all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%系統離散化及初始值設置%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fs=1000;%采樣頻率
N1=1000;%采樣點數
N=1000;
dfs=Fs/N;%頻率分辨率
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個實際頻率點
%sys=tf([-5],[2.0e-5 2e-9 1]);
sys=tf([3.355e7],[1 1.504e3 5.394e5 3.291e7]);
sysd=c2d(sys,1/Fs,'tustin');
[num,den]=tfdata(sysd,'v');
[h0,f0]=freqz(num,den,N/2,Fs);
mag=abs(h0);
ph=angle(h0);
ph=ph*180/pi;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%輸入信號%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=0.5;%頻率F1信號的幅度
A0=10;
f0=1;%起始頻率(Hz)
df=1;%頻率間隔
f1=30;%結束頻率
S=0;
x1=8*chirp(t1,1,N/Fs,60);%chirp
for i=1:N
x2(i)=8*sin(2*pi*(f0+(f1-f0)*i/N)*i/Fs);
end;
for i=N1:N+1 %補零消除柵欄效應
    S(i)=0;
    x(i)=0;
end;
[Ax,Px,A1x,P1x,A0x,P0x,y1]=qiuzhi(x1,N,num,den);
[Axf,Pxf,A1xf,P1xf,A0xf,P0xf,y1f]=qiuzhi(x2,N,num,den);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%基本繪圖%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
plot(t,x1,'r');
legend('組合正弦波');
grid;
title('輸入信號');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%FFT繪圖%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
plot(F,A1x,'r',F,A1xf,'b');
legend('chirp');
grid;%顯示換算后的FFT模值結果
title('輸入信號幅度-頻率曲線圖');
axis([0 200 -2 2]);