%fm.m
clear;
close all
clc;
t0=0.2; %信号的持续时间,用来定义时间向量
ts=0.001; %抽样间隔
fs=1/ts; %抽样频率
%*********************************************************************
fc=300; %载波频率,fc可以任意改变
%*********************************************************************
t=[-t0/2:ts:t0/2]; %时间向量
kf=100; %偏差常数
df=0.25; %所需的频率分辨率,用在求傅里叶变换时,它表示FFT的最小频率间隔
%*********************************************************************
m=sinc(100*t);
%m=sin(100*t); %调制信号,m(t)可以更改
%*********************************************************************
%*********************************************************************
%求信号m(t)的积分,对离散信号求积分采用求和的方式
int_m(1)=0;
for i=1:length(t)-1
int_m(i+1)=int_m(i)+m(i)*ts;
end
%*********************************************************************
%求傅里叶变换和计算已调信号
n1=length(m); %缩放,便于在频谱图上整体观察
M=fft(m,n1); %对调制信号m(t)求傅里叶变换
M=M/fs;
f1=[0:df:df*(length(m)-1)]-fs/2;
u=cos(2*pi*fc*t+2*pi*kf*int_m); %调制后的信号
f2=[0:df:df*(length(u)-1)]-fs/2;
n2=length(u);
U=fft(u,n2); %对调制后的信号u求傅里叶变换
U=U/fs; %缩放
%********************************************************************
%*********************************************************************
%经过高斯白噪声
%第一种加噪声方法,利用经典算法
snr_in_dB=25;
k=1;Eb=1e-3;
snr=10^(snr_in_dB/10);
sgma_idd=sqrt((Eb/snr)*fs/2);
awgn_noise=zeros(1,length(u)) ;
sign_noise=zeros(1,length(u)) ;
for i=1:length(u)
awgn_noise(i)=gngauss(sgma_idd);
end
u_add=u+awgn_noise;
%*********************************************************************
%第二种加噪声方法,直接利用matlab函数
%u_add=awgn(u,15);
%*********************************************************************
%不加噪声
%u_add=u;
%*********************************************************************
%*********************************************************************
%解调过程,通过希尔伯特变换求出正交部分,然后与载波抵消,得到相位
t1=[0:ts:ts*(length(u_add)-1)];
z=hilbert(u_add);
yq =z.*exp(-j*2*pi*fc*t1);
dem_m(1)=0;
dem_m(1,2:length(u_add))=diff(unwrap(angle(yq)));
dem =(1/(2*pi*kf))*dem_m;
%*********************************************************************
%*********************************************************************
%时域图,频谱图等
figure(1);
subplot(2,1,1);
plot(t,u(1:length(t)));
xlabel('时间t');
title('已调信号的时域图');
subplot(2,1,2);
plot(t,u_add(1:length(t)));
xlabel('时间t');
title('加高斯白噪声的已调信号时域图');
disp('按任意键可以看到原调制信号和已调信号的曲线');
pause
figure(2);
subplot(2,1,1);
plot(t,m(1:length(t)));
xlabel('时间t');
title('原调制信号的时域图');
subplot(2,1,2);
plot(t,u(1:length(t)));
xlabel('时间t');
title('已调信号的时域图');
disp('按任意键可以看到原调制信号和已调信号的在频域内的图形');
pause
figure(3);
subplot(2,1,1);
plot(f1,abs(fftshift(M)));%fftshift:将直流分量移到频谱中心
xlabel('频率f')
title('原调制信号的频谱图');
subplot(2,1,2);
plot(f2,abs(fftshift(U)));
xlabel('频率f');
title('已调信号的频谱图');
disp('按任意键可以看到在无噪声情况下的原调制信号的波形和解调信号的输出情况');
pause
figure(4);
subplot(2,1,1);
plot(t,m(1:length(t)));
xlabel('时间t');
title('原调制信号的时域图');
subplot(2,1,2);
plot(t,dem(1:length(t)));
xlabel('时间t');
title('解调后信号的时域波形');
评论1