echo?off
close?all
clear?all
clc

%%??FM調制
ps?=?625000;????????????????????%?數據率625KHz
Fs?=?80*10^6;???????????????????%?80M采樣率
dt?=?1/Fs;??????????????????????%?采樣間隔
a??=?1000;??????????????????????
t??=?0:dt:（a*Fs/ps-1）/Fs;???????%?產生時間序列
N??=?length（t）-1;
am?=?5;?????????????????????????%?設定調制信號幅度
fm?=?15*10^2;???????????????????%?調制信號頻率
fc?=?5*10^6;????????????????????%?載波頻率5M
mt?=?am*sin（2*pi*fm*t）;?????????%?生成調制信號
ct?=?cos（2*pi*fc*t）;????????????%?生成載波
kf?=?10000;?????????????????????%?調頻系數??頻偏75kHz

integral_mt（1）=0;
for?i=1:N???????????????????????%積分
????integral_mt（i+1）=integral_mt（i）+mt（i）*dt;
end
figure（100）;plot（tintegral_mt）;?title（‘integral_mt‘）;
sfm=?am*cos（2*pi*fc*t+2*pi*kf*integral_mt）;?%已調信號

%%%%%%%%%%%%畫圖%%%%%%%%%%%%
figure（1）;
subplot（311）;
plot（tmt）;?title（‘調制信號時域‘）;
axis（[0?10^（-3）?-5?5]）;??
subplot（312）;plot（tct）;?title（‘載波信號時域‘）;
axis（[0?10^（-5）?-1?2]）;??
subplot（313）;plot（tsfm）;title（‘已調信號時域‘）;
axis（[0?10^（-5）?-5?5]）;??

%%%%%%%%%%%%畫圖%%%%%%%%%%%%
figure（2）;
subplot（311）;
plot（（N/2:N/2）*Fs/N20*log10（abs（fftshift（fft（mt）））））;
title（‘調制信號頻域‘）;?
subplot（312）;
plot（（N/2:N/2）*Fs/N20*log10（abs（fftshift（fft（ct）））））;
title（‘載波信號頻譜‘）;?
subplot（313）;
plot（（N/2:N/2）*Fs/N20*log10（abs（fftshift（fft（sfm）））））;
title（‘已調信號頻域‘）;

%%%%%%%%%%%%%%%%%%生成IQ信號%%%%%%%%%%%%%%%%%%%%%
st=imag（hilbert（ct））;??%希爾伯特變化
ct=real（hilbert（ct））;??%cos（wc*t）
for?i=1:length（sfm）
Ih（i）=sfm（i）*ct（i）;
Qh（i）=sfm（i）*st（i）;
end

%%%%%%%%%%%%%%%%%%%%%濾波降頻%%%%%%%%%%%%%%%%%%%
Flp=fc;
b=fir1（200Flp/Fs‘low‘）;?%設計濾波器
%?figure（4）;??
%?stem（b）;????????????????%濾波器單位脈沖響應
Qm=filter（b1Qh）;?
Im=filter（b1Ih）;
%?figure（6）
%?subplot（211）;
plot（（N/2:N/2）*Fs/Nabs（fftshift（fft（Qm））））;
title（‘Q信號波形頻域‘）;?
%?subplot（212）;
plot（（N/2:N/2）*Fs/Nabs（fftshift（fft（Im））））;
title（‘I信號波形頻域‘）;?
%?figure（5）
%?subplot（211）;
plot（（N/2:N/2）*Fs/Nabs（ff