最新AM仿真程序.docx
《最新AM仿真程序.docx》由会员分享,可在线阅读,更多相关《最新AM仿真程序.docx(15页珍藏版)》请在冰点文库上搜索。
最新AM仿真程序
AM仿真程序
仿真程序
PM
首先任意给定一个已知调制信号m(t)=cos(pi*10*t),进行相位调制时要用到傅里叶变换,因此先编写傅里叶变换的m文件用作主函数调用,其m文件代码如下:
%求傅里叶变换的子函数
function[M,m,df]=fftseq(m,ts,df)
fs=1/ts;
ifnargin==2n1=0;%nargin为输入参量的个数
elsen1=fs/df;
end
n2=length(m);
n=2^(max(nextpow2(n1),nextpow2(n2)));
%nextpow2(n)取n最接近的较大2次幂
M=fft(m,n);
%M为信号m的傅里叶变换,n为快速傅里叶变换的点数,及基n-FFT变换
m=[m,zeros(1,n-n2)];%构建新的m信号
df=fs/n;%重新定义频率分辨率
上述m文件以“fftseq.m”保存。
在实现相位解调时要调用两个子函数,分述如下:
%求信号相角的子函数,这是调频、调相都要用到的方法
function[v,phi]=env_phas(x,ts,f0)
ifnargout==2%nargout为输出变数的个数
z=loweq(x,ts,f0);%产生调制信号的正交分量
phi=angle(z);%angle是对一个复数求相角的函数
end
v=abs(hilbert(x));%abs用来求复数hilbert(x)的模
上述m文件以“env_phas.m”保存。
%产生调制信号的正交分量
functionx1=loweq(x,ts,f0)
t=[0:
ts:
ts*(length(x)-1)];
z=hilbert(x);%希尔伯特变换对的利用---通过实部来求虚部
x1=z.*exp(-j*2*pi*f0*t);
%产生信号z的正交分量,%并将z信号与它的正交分量加在一起
上述m文件以“loweq.m”保存
%主程序
t0=1;%信号的持续时间,用来定义时间向量
ts=0.001;%抽样间隔
fs=1/ts;%抽样频率
fc=100;%载波频率,fc可以任意改变
t=[-t0/2:
ts:
t0/2];%时间向量
kf=100;%偏差常数
df=0.25;
%所需的频率分辨率,用在求傅里叶变换时,它表示FFT的最小频率间隔
m=cos(pi*10*t);%调制信号,m(t)可以任意更改
int_m
(1)=0;%求信号m(t)的积分
fori=1:
length(t)-1
int_m(i+1)=int_m(i)+m(i)*ts;
end
[M,m,df1]=fftseq(m,ts,df);%对调制信号m(t)求傅里叶变换
M=M/fs;%缩放,便于在频谱图上整体观察
f=[0:
df1:
df1*(length(m)-1)]-fs/2;%时间向量对应的频率向量
u=cos(2*pi*fc*t+2*pi*kf*int_m);%调制后的信号
[U,u,df1]=fftseq(u,ts,df);%对调制后的信号u求傅里叶变换
U=U/fs;%缩放
%通过调用子程序env_phas和loweq来实现解调功能
[v,phase]=env_phas(u,ts,fc);%解调,求出u的相位
phi=unwrap(phase);
%校正相位角,使相位在整体上连续,便于后面对该相位角求导dem=(1/(2*pi*kf))*(diff(phi)*fs);%对校正后的相位求导%再经一些线性变换来恢复原调制信号%乘以fs是为了恢复原信号,因为前面使用了缩放
subplot(2,3,1)%子图形式显示结果
plot(t,m(1:
length(t)))%现在的m信号是重新构建的信号,%因为在对m求傅里叶变换时m=[m,zeros(1,n-n2)]
axis([-0.50.5-11])%定义两轴的刻度
xlabel('时间t')
title('原调制信号的时域图')
subplot(2,3,4)
plot(t,u(1:
length(t)))
axis([-0.50.5-11])
xlabel('时间t')
title('已调信号的时域图')
subplot(2,3,2)
plot(f,abs(fftshift(M)))%fftshift:
将FFT中的DC分量移到频谱中心
axis([-60060000.05])
xlabel('频率f')
title('原调制信号的频谱图')
subplot(2,3,5)
plot(f,abs(fftshift(U)))
axis([-60060000.05])
xlabel('频率f')
title('已调信号的频谱图')
subplot(2,3,3)
plot(t,m(1:
length(t)))
axis([-0.50.5-11])
xlabel('时间t')
title('原调制信号的时域图')
subplot(2,3,6)
plot(t,m(1:
length(t)))
axis([-0.50.5-11])
xlabel('时间t')
title('解调后信号的时域波形')
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;
fori=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解调*******************
fori=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;
ifnargin==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;
ifnargin==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;
fori=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;%生成含高斯白躁声的已调信号(信号通
%过信道传输)
fori=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;%生成含高斯白躁声的已调信号(信号通过信道传输)
fori=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('含大信噪比高斯白噪声解调信号的时域图');
%******************结束*******************
AM
fm=100;
fc=500;
fs=5000;
Am=1;
A=2;
N=512;
K=N-1;
n=0:
N-1;
t=(0:
1/fs:
K/fs);
yt=Am*cos(2*pi*fm*t);
figure
(1)
subplot(1,1,1),plot(t,yt),title('频率为3000的调制信号f1的时时域波');
y0=A+yt;
y2=y0.*cos(2*pi*fc*n/fs);
y3=fft(y2,N);%fft变换
q1=(0:
N/2-1)*fs/N;
mx1=abs(y3(1:
N/2));
figure
(2)
subplot(2,1,1);
plot(t,y2);
title('已调信号的时时域波');
subplot(2,1,2);
plot(q1,mx1);
title('f1已调信号的频谱');%绘图
yc=cos(2*pi*fc*t);
figure(3)
subplot(2,1,1),plot(t,yc),title('载波fc时域波形')
N=512;
n=0:
N-1;
yc1=Am*cos(2*pi*fc*n/fs);
y3=fft(yc1,N);
q=(0:
N/2-1)*fs/N;
mx=abs(y3(1:
N/2));
figure(3)
subplot(2,1,2),plot(q,mx),title('载波fc频谱')
y4=0.01*randn(1,length(t));%用RANDN产生高斯分布序列
w=y4.^2;%噪声功率
figure(4)
subplot(2,1,1);
plot(t,y4);
title('高斯白噪声时域波形')
y5=fft(y4,N);
q2=(0:
N/2-1)*fs/N;
mx2=abs(y5(1:
N/2));
subplot(2,1,2),plot(q2,mx2),title('高斯白噪声频域波形')
y6=y2+y4;
figure(5)
subplot(2,1,1),plot(t,y6),title('叠加后的调制信号时域波形')
q3=q1;
mx3=mx1+mx2;
subplot(2,1,2),plot(q3,mx3),title('叠加后的调制信号频谱波形')%调制
yv=y6.*yc;%乘以载波进行解调
Ws=yv.^2;
p1=fc-fm;
[k,Wn,beta,ftype]=kaiserord([p1fc],[10],[0.050.01],fs);%Fir数字低通滤波
window=kaiser(k+1,beta);%使用kaiser窗函数
b=fir1(k,Wn,ftype,window,'noscale');%使用标准频率响应的加窗设计函数
yt=filter(b,1,yv);
yssdb=yt.*2-2;
figure(6)
subplot(2,1,1),plot(t,yssdb),title('经过低通已调信号的时域波形采样')
y9=fft(yssdb,N);
q=(0:
N/2-1)*fs/N;
mx=abs(y9(1:
N/2));
subplot(2,1,2),plot(q,mx),title('经过低通已调信号频域波形')%解调
ro=y9-yt;
W=(yt.^2).*(1/2);
R=W/w
r=W/ro
G=r/R