基于Matlab 的FM 调制解调仿真
叶傻逼 白痴NO.1
1.1 FM 调制模型的建立
图2 FM 调制模型
其中，()m t 为基带调制信号，设调制信号为
()cos(2)m m t A f t π=
设正弦载波为
()cos(2)c c t f t π=
信号传输信道为高斯白噪声信道，其功率为2σ。
1.2 调制过程分析
在调制时，调制信号的频率去控制载波的频率的变化，载波的瞬时频偏随调制信号()m t 成正比例变化，即
(
)
()f d t K m t dt
?=
式中，f K 为调频灵敏度（()rad s V ?）。 这时相位偏移为
()()f t K m d ?ττ=?
则可得到调频信号为
()cos ()FM c f s t A t K m d ωττ??=+??
?
调制信号产生的M 文件：
dt=0.001; %设定时间步长
t=0:dt:1.5; %产生时间向量
am=15; %设定调制信号幅度←可更改 fm=15; %设定调制信号频率←可更改 mt=am*cos(2*pi*fm*t); %生成调制信号
fc=50; %设定载波频率←可更改 ct=cos(2*pi*fc*t); %生成载波 kf=10; %设定调频指数 int_mt(1)=0; %对mt 进行积分 for i=1:length(t)-1
int_mt(i+1)=int_mt(i)+mt(i)*dt; end
sfm=am*cos(2*pi*fc*t+2*pi*kf*int_mt); %调制，产生已调信号
0.5
1 1.5
-10
010时间t
调制信号的时域图
0.5
1
1.5
-101时间t
载波的时域图
00.5
1 1.5
-10
010时间t
已调信号的时域图
图3 FM 调制
1.3 FM 解调模型的建立
调制信号的解调分为相干解调和非相干解调两种。相干解调仅仅适用于窄带调频信号，且需同步信号，故应用范围受限；而非相干解调不需同步信号，且对于NBFM 信号和WBFM 信号均适用，因此是FM 系统的主要解调方式。在本仿真的过程中我们选择用非相干解调方法进行解调。
图4 FM 解调模型
非相干解调器由限幅器、鉴频器和低通滤波器等组成，其方框图如图5所示。限幅器输入为已调频信号和噪声，限幅器是为了消除接收信号在幅度上可能出现的畸变；带通滤波器的作用是用来限制带外噪声，使调频信号顺利通过。鉴频器中的微分器把调频信号变成调幅调频波，然后由包络检波器检出包络，最后通过低通滤波器取出调制信号。
1.4 解调过程分析
设输入调频信号为
()()cos(())t
t FM c f S t S t A t K m d ωττ-∞
==+?
微分器的作用是把调频信号变成调幅调频波。微分器输出为
[]()()()()sin(())
i FM d t
c f c f dS t dS t S t dt dt
K m t t K m d ωωττ+-∞==
=-+?
包络检波的作用是从输出信号的幅度变化中检出调制信号。包络检波器输出为
[]()()()o d c f d c d f S t K Km t K K Km t ωω++==
d K 称为鉴频灵敏度（V Hz ）
，是已调信号单位频偏对应的调制信号的幅度，经低通滤波器后加隔直流电容，隔除无用的直流，得
()()o d f m t K K m t =
微分器通过程序实现，代码如下：
for i=1:length(t)-1 %接受信号通过微分器处理 diff_nsfm(i)=(nsfm(i+1)-nsfm(i))./dt; end
diff_nsfmn = abs(hilbert(diff_nsfm)); %hilbert 变换，求绝对值得到瞬时幅度
（包络检波）
通过M 文件绘制出两种不同信噪比解调的输出波形如下：
00.51 1.5
-5
5时间t
调制信号的时域图
0.5
1
1.5
-50
5时间t
无噪声条件下已调信号的时域图
00.5
1 1.5
-5
5时间t
无噪声条件下解调信号的时域图
图5 FM 解调
1.5 高斯白噪声信道特性
设正弦波通过加性高斯白噪声信道后的信号为
()cos()()c r t A t n t ωθ=++
其中，白噪声()n t 的取值的概率分布服从高斯分布。
MATLAB 本身自带了标准高斯分布的内部函数randn 。randn 函数产生的随机序列服从均值为0m =，方差21σ=的高斯分布。
正弦波通过加性高斯白噪声信道后的信号为
()cos()()c r t A t n t ωθ=++
故其有用信号功率为
2
2
A S =
噪声功率为
2N σ=
信噪比S
N
满足公式
1010log ()S B N
=
则可得到公式
22
10
210
B A σ=
?
我们可以通过这个公式方便的设置高斯白噪声的方差。
在本仿真过程中，我们选择了10db 和30db 两种不同信噪比以示区别，其时域图如图7和图8。
00.51 1.5
-5
5时间t
调制信号的时域图
0.5
1
1.5
-50
5时间t
无噪声条件下已调信号的时域图
00.5
1 1.5
-5
5时间t
无噪声条件下解调信号的时域图
图6 无噪声条件下已调信号的时域图
00.51 1.5
-5
5时间t
调制信号的时域图
500
1000
1500
-5000
5000
时间t
含小信噪比高斯白噪声已调信号的时域图
00.5
1 1.5
-20
20时间t
含小信噪比高斯白噪声解调信号的时域图
图7 含小信噪比高斯白噪声已调信号的时域图
00.51 1.5
-5
5时间t
调制信号的时域图
500
1000
1500
-5000
5000
时间t
含大信噪比高斯白噪声已调信号的时域图
00.5
1 1.5
-5
5时间t
含大信噪比高斯白噪声解调信号的时域图
图8 含大信噪比高斯白噪声已调信号的时域图
三．仿真实现
图10 程序流程图
3.1 MATLAB源代码
%FM调制解调系统.m
%频率调制与解调的Matlab演示源程序
%可以任意改原调制信号函数m(t)
%信息工程陈丽丹07323202 %·*·*·*·*·*·*·*·*·*·*·*·*·*·*·*·
%*****************初始化******************
echo off
close all
clear all
clc
%***************************************** %·*·*·*·*·*·*·*·*·*·*·*·*·*·*·*·
%****************FM调制*******************
dt=0.001; %设定时间步长
t=0:dt:1.5; %产生时间向量
am=5; %设定调制信号幅度
fm=5; %设定调制信号频率
mt=am*cos(2*pi*fm*t); %生成调制信号
fc=50; %设定载波频率
ct=cos(2*pi*fc*t); %生成载波
kf=10; %设定调频指数
int_mt(1)=0;
for i=1:length(t)-1
int_mt(i+1)=int_mt(i)+mt(i)*dt; %求信号m(t)的积分
end %调制，产生已调信号
sfm=am*cos(2*pi*fc*t+2*pi*kf*int_mt); %调制信号
%***************************************** %·*·*·*·*·*·*·*·*·*·*·*·*·*·*·*·
%*************添加高斯白噪声**************
sn1=10; %设定信躁比(小信噪比)
sn2=30; %设定信躁比(大信噪比)
sn=0; %设定信躁比(无信噪比)
db=am^2/(2*(10^(sn/10))); %计算对应的高斯白躁声的方差
n=sqrt(db)*randn(size(t)); %生成高斯白躁声
nsfm=n+sfm; %生成含高斯白躁声的已调信号（信号通
%过信道传输）
%***************************************** %·*·*·*·*·*·*·*·*·*·*·*·*·*·*·*·
%****************FM解调*******************
for i=1:length(t)-1 %接受信号通过微分器处理diff_nsfm(i)=(nsfm(i+1)-nsfm(i))./dt;
end
diff_nsfmn = abs(hilbert(diff_nsfm)); %hilbert变换，求绝对值得到瞬时幅度（包络检波）zero=(max(diff_nsfmn)-min(diff_nsfmn))/2;
diff_nsfmn1=diff_nsfmn-zero;
%***************************************** %·*·*·*·*·*·*·*·*·*·*·*·*·*·*·*·
%**************时域到频域转换**************
ts=0.001; %抽样间隔
fs=1/ts; %抽样频率
df=0.25; %所需的频率分辨率，用在求傅里叶变换
%时，它表示FFT的最小频率间隔
%*****对调制信号m(t)求傅里叶变换*****
m=am*cos(2*pi*fm*t); %原调信号
fs=1/ts;
if nargin==2
n1=0;
else
n1=fs/df;
end
n2=length(m);
n=2^(max(nextpow2(n1),nextpow2(n2)));
M=fft(m,n);
m=[m,zeros(1,n-n2)];
df1=fs/n; %以上程序是对调制后的信号u求傅里变换
M=M/fs; %缩放，便于在频铺图上整体观察
f=[0:df1:df1*(length(m)-1)]-fs/2; %时间向量对应的频率向量
%************对已调信号u求傅里变换**********
fs=1/ts;
if nargin==2
n1=0;
else
n1=fs/df;
end
n2=length(sfm);
n=2^(max(nextpow2(n1),nextpow2(n2)));
U=fft(sfm,n);
u=[sfm,zeros(1,n-n2)];
df1=fs/n; %以上是对已调信号u求傅里变换
U=U/fs; %缩放
%******************************************
%***************************************** %·*·*·*·*·*·*·*·*·*·*·*·*·*·*·*·
%***************显示程序******************
disp('按任意键可以看到原调制信号、载波信号和已调信号的曲线')
pause
%**************figure(1)******************
figure(1)
subplot(3,1,1);plot(t,mt); %绘制调制信号的时域图
xlabel('时间t');
title('调制信号的时域图');
subplot(3,1,2);plot(t,ct); %绘制载波的时域图
xlabel('时间t');
title('载波的时域图');
subplot(3,1,3);
plot(t,sfm); %绘制已调信号的时域图
xlabel('时间t');
title('已调信号的时域图');
%******************************************
disp('按任意键可以看到原调制信号和已调信号在频域内的图形')
pause
%************figure(2)*********************
figure(2)
subplot(2,1,1)
plot(f,abs(fftshift(M))) %fftshift:将FFT中的DC分量移到频谱中心xlabel('频率f')
title('原调制信号的频谱图')
subplot(2,1,2)
plot(f,abs(fftshift(U)))
xlabel('频率f')
title('已调信号的频谱图')
%******************************************
disp('按任意键可以看到原调制信号、无噪声条件下已调信号和解调信号的曲线')
pause
%**************figure(3)******************
figure(3)
subplot(3,1,1);plot(t,mt); %绘制调制信号的时域图
xlabel('时间t');
title('调制信号的时域图');
subplot(3,1,2);plot(t,sfm); %绘制已调信号的时域图
xlabel('时间t');
title('无噪声条件下已调信号的时域图');
nsfm=sfm;
for i=1:length(t)-1 %接受信号通过微分器处理diff_nsfm(i)=(nsfm(i+1)-nsfm(i))./dt;
end
diff_nsfmn = abs(hilbert(diff_nsfm)); %hilbert变换，求绝对值得到瞬时幅度（包络检波）zero=(max(diff_nsfmn)-min(diff_nsfmn))/2;
diff_nsfmn1=diff_nsfmn-zero;
subplot(3,1,3); %绘制无噪声条件下解调信号的时域图
plot((1:length(diff_nsfmn1))./1000,diff_nsfmn1./400,'r');
xlabel('时间t');
title('无噪声条件下解调信号的时域图');
%*****************************************
disp('按任意键可以看到原调制信号、小信噪比高斯白噪声条件下已调信号和解调信号已调信号的曲线')
pause
%**************figure(4)******************
figure(4)
subplot(3,1,1);plot(t,mt); %绘制调制信号的时域图
xlabel('时间t');
title('调制信号的时域图');
db1=am^2/(2*(10^(sn1/10))); %计算对应的小信噪比高斯白躁声的方差n1=sqrt(db1)*randn(size(t)); %生成高斯白躁声
nsfm1=n1+sfm; %生成含高斯白躁声的已调信号（信号通
%过信道传输）
for i=1:length(t)-1 %接受信号通过微分器处理diff_nsfm1(i)=(nsfm1(i+1)-nsfm1(i))./dt;
end
diff_nsfmn1 = abs(hilbert(diff_nsfm1)); %hilbert变换，求绝对值得到瞬时幅度（包络检波）zero=(max(diff_nsfmn)-min(diff_nsfmn))/2;
diff_nsfmn1=diff_nsfmn1-zero;
subplot(3,1,2);
plot(1:length(diff_nsfm),diff_nsfm); %绘制含小信噪比高斯白噪声已调信号的时域图xlabel('时间t');
title('含小信噪比高斯白噪声已调信号的时域图');
subplot(3,1,3); %绘制含小信噪比高斯白噪声解调信号的时域图plot((1:length(diff_nsfmn1))./1000,diff_nsfmn1./400,'r');
xlabel('时间t');
title('含小信噪比高斯白噪声解调信号的时域图');
%*****************************************
disp('按任意键可以看到原调制信号、大信噪比高斯白噪声条件下已调信号和解调信号已调信号的曲线')
pause
%**************figure(5)******************
figure(5)
subplot(3,1,1);plot(t,mt); %绘制调制信号的时域图
xlabel('时间t');
title('调制信号的时域图');
db1=am^2/(2*(10^(sn2/10))); %计算对应的大信噪比高斯白躁声的方差n1=sqrt(db1)*randn(size(t)); %生成高斯白躁声
nsfm1=n1+sfm; %生成含高斯白躁声的已调信号（信号通过信道传输）
for i=1:length(t)-1 %接受信号通过微分器处理diff_nsfm1(i)=(nsfm1(i+1)-nsfm1(i))./dt;
end
diff_nsfmn1 = abs(hilbert(diff_nsfm1)); %hilbert变换，求绝对值得到瞬时幅度（包
%络检波）
zero=(max(diff_nsfmn)-min(diff_nsfmn))/2;
diff_nsfmn1=diff_nsfmn1-zero;
subplot(3,1,2);
plot(1:length(diff_nsfm1),diff_nsfm1); %绘制含大信噪比高斯白噪声已调信号
%的时域图
xlabel('时间t');
title('含大信噪比高斯白噪声已调信号的时域图');
subplot(3,1,3); %绘制含大信噪比高斯白噪声解调信号
%的时域图
plot((1:length(diff_nsfmn1))./1000,diff_nsfmn1./400,'r');
xlabel('时间t');
title('含大信噪比高斯白噪声解调信号的时域图');
%***************************************** %******************结束*******************
3.2 仿真结果