《AM仿真程序.doc》由会员分享，可在线阅读，更多相关《AM仿真程序.doc（14页珍藏版）》请在三一文库上搜索。
1、仿真程序PM首先任意给定一个已知调制信号m(t)=cos(pi*10*t), 进行相位调制时要用到傅里叶变换，因此先编写傅里叶变换的m文件用作主函数调用，其m文件代码如下：%求傅里叶变换的子函数function M,m,df=fftseq(m,ts,df)fs=1/ts;if nargin=2 n1=0; %nargin为输入参量的个数else n1=fs/df;endn2=length(m);n=2(max(nextpow2(n1),nextpow2(n2); %nextpow2(n)取n最接近的较大2次幂M=fft(m,n); %M为信号m的傅里叶变换，n为快速傅里叶变换的点数，及基n-F
2、FT变换m=m,zeros(1,n-n2); %构建新的m信号df=fs/n; %重新定义频率分辨率上述m文件以“fftseq.m”保存。 在实现相位解调时要调用两个子函数，分述如下：%求信号相角的子函数，这是调频、调相都要用到的方法function v,phi=env_phas(x,ts,f0)if nargout=2 %nargout为输出变数的个数 z=loweq(x,ts,f0); %产生调制信号的正交分量 phi=angle(z); %angle是对一个复数求相角的函数endv=abs(hilbert(x); %abs用来求复数hilbert(x)的模上述m文件以“env_phas.
3、m”保存。%产生调制信号的正交分量function x1=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; %所需的频率
4、分辨率，用在求傅里叶变换时，它表示FFT的最小频率间隔m=cos(pi*10*t); %调制信号，m(t)可以任意更改int_m(1)=0; %求信号m(t)的积分for i=1:length(t)-1 int_m(i+1)=int_m(i)+m(i)*ts;endM,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,
5、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求傅
6、里叶变换时m=m,zeros(1,n-n2)axis(-0.5 0.5 -1 1) %定义两轴的刻度xlabel(时间t) title(原调制信号的时域图)subplot(2,3,4)plot(t,u(1:length(t)axis(-0.5 0.5 -1 1)xlabel(时间t)title(已调信号的时域图)subplot(2,3,2)plot(f,abs(fftshift(M) %fftshift：将FFT中的DC分量移到频谱中心axis(-600 600 0 0.05)xlabel(频率f)title(原调制信号的频谱图)subplot(2,3,5)plot(f,abs(fftshif
7、t(U)axis(-600 600 0 0.05)xlabel(频率f)title(已调信号的频谱图)subplot(2,3,3)plot(t,m(1:length(t)axis(-0.5 0.5 -1 1)xlabel(时间t)title(原调制信号的时域图)subplot(2,3,6)plot(t,m(1:length(t)axis(-0.5 0.5 -1 1)xlabel(时间t)title(解调后信号的时域波形)FMdt=0.001; %设定时间步长t=0:dt:1.5; %产生时间向量am=5; %设定调制信号幅度fm=5; %设定调制信号频率mt=am*cos(2*pi*fm*t)
8、; %生成调制信号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=am2/(2*(10(sn/10); %计算对应的高斯白躁声的方差n
9、=sqrt(db)*randn(size(t); %生成高斯白躁声nsfm=n+sfm; %生成含高斯白躁声的已调信号（信号通%过信道传输）%*%*FM解调*for i=1:length(t)-1 %接受信号通过微分器处理 diff_nsfm(i)=(nsfm(i+1)-nsfm(i)./dt;enddiff_nsfmn = abs(hilbert(diff_nsfm); %hilbert变换，求绝对值得到瞬时幅度（包络检波）zero=(max(diff_nsfmn)-min(diff_nsfmn)/2;diff_nsfmn1=diff_nsfmn-zero;%*%*时域到频域转换*ts=0.
10、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;endn2=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
11、(m)-1)-fs/2; %时间向量对应的频率向量%*对已调信号u求傅里变换*fs=1/ts;if nargin=2 n1=0;else n1=fs/df;endn2=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); %绘制调制信号的时域图x
12、label(时间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(原调制信号
13、的频谱图)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)-
14、1 %接受信号通过微分器处理 diff_nsfm(i)=(nsfm(i+1)-nsfm(i)./dt;enddiff_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(无噪声条件下解调信
15、号的时域图);%*disp(按任意键可以看到原调制信号、小信噪比高斯白噪声条件下已调信号和解调信号已调信号的曲线)pause%*figure(4)*figure(4)subplot(3,1,1);plot(t,mt); %绘制调制信号的时域图xlabel(时间t);title(调制信号的时域图);db1=am2/(2*(10(sn1/10); %计算对应的小信噪比高斯白躁声的方差n1=sqrt(db1)*randn(size(t); %生成高斯白躁声nsfm1=n1+sfm; %生成含高斯白躁声的已调信号（信号通%过信道传输）for i=1:length(t)-1 %接受信号通过微分器处理 d
16、iff_nsfm1(i)=(nsfm1(i+1)-nsfm1(i)./dt;enddiff_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
17、); %绘制含小信噪比高斯白噪声解调信号的时域图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=am2/(2*(10(sn2/10); %计算对应的大信噪比高斯白躁声的方差
18、n1=sqrt(db1)*randn(size(t); %生成高斯白躁声nsfm1=n1+sfm; %生成含高斯白躁声的已调信号（信号通过信道传输）for i=1:length(t)-1 %接受信号通过微分器处理 diff_nsfm1(i)=(nsfm1(i+1)-nsfm1(i)./dt;enddiff_nsfmn1 = abs(hilbert(diff_nsfm1); %hilbert变换，求绝对值得到瞬时幅度（包%络检波）zero=(max(diff_nsfmn)-min(diff_nsfmn)/2;diff_nsfmn1=diff_nsfmn1-zero;subplot(3,1,2);
19、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(含大信噪比高斯白噪声解调信号的时域图);%*结 束*AMfm=100;fc=500;fs=5000;Am=1;A=2;N=512;K=N-1;n=0:N-1;t=(0:1/fs:K/
20、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);fi
21、gure(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
22、);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(p1 fc,1 0,0.05 0.01,fs)
23、; %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