close all
t0=8;  %信號的持續時間,用來定義時間軸的取值范圍
ts=0.0001;  %抽樣間隔
fs=1/ts;  %抽樣頻率
fc=10000;  %載波頻率
t=[-t0/2:ts:t0/2];  %時間軸的取值
kf=5;   %壓控振蕩系數
df=0.25;  %所需的頻率分辨率，用在求傅里葉變換時，它表示FFT的最小頻率間隔
m=2*pi*t;  %調制信號
%*********************************************************************
jifen(1)=0;  %求信號m(t)的積分
for i=1:length(t)-1  
    jifen(i+1)=jifen(i)+m(i)*ts;
end
%*********************************************************************
[M,m,df1]=fourier(m,ts,df); %對調制信號m(t)求傅里葉變換，需調用子函數fourier
M=M/fs;    %縮放，便于在頻譜圖上整體觀察
f=[1:length(m)]-fs/2;  %時間軸上對應的頻率向量
u=1.414*cos(2*pi*fc*t+2*pi*kf*jifen);    %調制后的信號
[U,u,df1]=fourier(u,ts,df);       %對調制后的信號u(t)求傅里葉變換
U=U/fs;    %縮放
%*********************************************************************
t1=[0:ts:ts*(length(u)-1)];  %對已調信號進行解調
z=hilbert(u);  %希爾伯特變換
yq =z.*exp(-j*2*pi*fc*t1);
dem =diff(unwrap(angle(yq))*fs)/(2*pi*kf);
                 %unwrap函數用于恢復相位并去掉2π相位疊加的影響
for i=1:(length(t)-1)/2  
    dem(i)=-dem(i);
end
%*********************************************************************
figure(1)
subplot(2,1,1)
plot(t,m(1:length(t))) %此時的m信號是重新構建的信號，
                       %因為在對m求傅里葉變換時m=[m,zeros(1,n-n2)]
xlabel('時間t')
title('調制信號的時域波形')
subplot(2,1,2)
plot(t,u(1:length(t)))
axis([-1.5,1.5,-1.5,1.5]);
xlabel('時間t')
title('已調信號的時域波形')

figure(2)
subplot(2,1,1)
plot(f,abs(fftshift(M))) %fftshift：將FFT中的DC分量移到頻譜中心
axis([58000,63000,0,2]);
xlabel('頻率f')
title('調制信號的頻譜圖')
subplot(2,1,2)
plot(f,abs(fftshift(U)))
axis([58000,63000,0,0.5]);
xlabel('頻率f')
title('已調信號的頻譜圖')

figure(3)
subplot(2,1,1)
plot(t,m(1:length(t)))
xlabel('時間t')
title('調制信號的時域波形')
subplot(2,1,2)
plot(t,dem(1:length(t)))
axis([-t0/2,t0/2,-t0*5,t0*5]);
xlabel('時間t')
title('解調信號的時域波形')