北邮通信原理编程实验完整代码.doc

北邮通信原理编程实验完整代码告各位学弟学妹书：此实验对于大多数人来说甚是烦人，但是我还是要劝你：尽力去学吧，学的时候觉得没什么用，那是因为不了解其博大精深，如果你一不小心顶住压力认真学了（譬如me），日后涉及相关领域时，你就能够触类旁通，瞬间迸发出灵感，撞击出智慧的火花，这种感觉就好像搂着一群漂亮的姑娘行走在春风中一样，原来那些朦胧的理论都变成了清晰的线条，串在一起：你会有一种类似于愚公般的给子子孙孙签订卖身契投身于学术的激进想法，不过这种极左的思想毕竟要不得除非你想女朋友和你分手，于激励学习便是极好的。实验一：/基带信号m(t)=sin(2000*pi*t)+2cos(1000*pi*t),fc=20khz,求AM,DSB-SC,SSBclear all exec t2f.sci; exec f2t.sci;fs=800;/采样速率T=200;/截短时间N=T*fs;/采样点数dt=1/fs; /时域采样间隔t=-T/2:dt:T/2-dt; /时域采样点df=1/T; /频域采样间隔f=-fs/2:df:fs/2-df; /频域采样点数fm1=1; /待观测正弦波频率，单位KHz，下同fm2=0.5; /待观测余弦波频率fc=20; /载波频率/以上为初始化参数设置m1=sin(2*%pi)*fm1*t); /待观测正弦波部分M1=t2f(m1,fs); /傅里叶变换MH1=-%i*sign(f).*M1; /希尔伯特变换mh1=real(f2t(MH1,fs); /希尔伯特反变换m2=2*cos(2*%pi)*fm2*t); /待观测余弦波部分M2=t2f(m2,fs); /傅里叶变换MH2=-%i*sign(f).*M2; /希尔伯特变换mh2=real(f2t(MH2,fs); /希尔伯特反变换s1=(1+(m1+m2)/abs(max(m1+m2).*cos(2*%pi)*fc*t); /AM信号时域表达式S1=t2f(s1,fs); /AM信号频域表达式s2=(m1+m2).*cos(2*%pi)*fc*t); /DSB-SC信号时域表达式S2=t2f(s2,fs); /DSB-SC信号频域表达式s3=(m1+m2).*cos(2*%pi)*fc*t)-(mh1+mh2).*sin(2*%pi)*fc*t); /SSB信号时域表达式，以上边带为例S3=t2f(s3,fs); /SSB信号上边带频域表达式/以上是仿真计算部分/以下为绘图部分/AM信号xset(window,1)plot(f,abs(S1)title(AM信号频谱)xlabel(f)ylabel(S(f)mtlb_axis(-25,25,0,max(abs(S1);xset(window,2)plot(t,s1)title(AM信号波形)xlabel(t)ylabel(s(t)mtlb_axis(-3,3,-3,3);/DSB-SC信号xset(window,3)plot(f,abs(S2)title(DSB-SC信号频谱)xlabel(f)ylabel(S(f)mtlb_axis(-25,25,0,max(abs(S2);xset(window,4)plot(t,s2)title(DSB-SC信号波形)xlabel(t)ylabel(s(t)mtlb_axis(-1,4,-3,3);/SSB信号（以上边带为例）xset(window,5)plot(f,abs(S3)title(SSB信号频谱)xlabel(f)ylabel(S(f)mtlb_axis(-25,25,0,max(abs(S3)xset(window,6)plot(t,s3)title(SSB信号波形)xlabel(t)ylabel(s(t)mtlb_axis(0,6,-3,3)实验二：/m(t)=sin(2000*pi*t)+2cos(1000*pi*t)+4sin(500*pi*t+pi/3),fc=40KHz,Kf=5kHz/V，仿真产生FM，观察频谱和波形/初始化参数设置clear all;exec t2f.sci;exec f2t.sci;xdel(winsid(); /关闭所有图形窗口fs=4000; /采样速率T=16; /截短时间N=T*fs;/采样点数dt=1/fs; /时域采样间隔t=-T/2:dt:T/2-dt; /时域采样点df=1/T; /频域采样间隔f=-fs/2:df:fs/2-df; /频域采样点数fm1=1;fm2=0.5;fm3=0.25;fc=40;Kf=5; /频偏常数/仿真计算m1=sin(2*%pi)*fm1*t);phi1=(2*%pi)*Kf*cumsum(m1)*dt; /cumsum函数为求累积和，相当于积分m2=2*cos(2*%pi)*fm2*t);phi2=(2*%pi)*Kf*cumsum(m2)*dt;m3=4*sin(2*%pi)*fm3*t+%pi/3);phi3=(2*%pi)*Kf*cumsum(m3)*dt;s=cos(2*%pi)*fc*t+phi1+phi2+phi3); /FM信号时域表达式S=t2f(s,fs); /FM信号频域表达式m=sin(2*%pi*fm1*t)+2*cos(2*%pi*fm2*t)+4*sin(2*%pi*fm3*t+%pi/3); /调制信号表达式M=t2f(m,fs);/绘图/FM已调信号xset(window,1)plot(f,abs(S)title(FM信号频谱)xlabel(f)ylabel(S(f)mtlb_axis(-80,80,0,max(abs(S);xset(window,2)plot(t,s)title(FM信号波形)xlabel(t)ylabel(s(t)mtlb_axis(0,3,-2,2);/调制信号xset(window,3)plot(f,abs(M)title(m(t)信号频谱)xlabel(f)ylabel(S(f)mtlb_axis(-5,5,0,max(abs(M);xset(window,4)plot(t,m)title(m(t)信号波形)xlabel(t)ylabel(m(t)mtlb_axis(-8,8,-6,8);实验三：/通过仿真测量占空比为50%、75%以及100%的单双极性归零码波形及其功率/谱，分析不同占空比对仿真结果的影响clear all;xdel(winsid(); /关闭所有图形窗口exec t2f.sci; exec f2t.sci;L=32; /每个码元间隔内的采样点数N=213; /总采样点数M=N/L; /总码元数Rb=2; /码元速率Ts=1/Rb; /比特间隔 fs=L/Ts; /采样速率T=N/fs; /截短时间 Bs=fs/2; /系统带宽 t=-T/2+0:N-1/fs; /时域采样点f=-Bs+0:N-1/T; /频域采样点L0=input(请输入占空比(01)：)EP=zeros(1,N);ch=input(请选择要观察的码型：1-单极性；2-双极性)for loop=1:1000 /1000次样本函数取平均 if ch=1 /单极性条件 a=(rand(1,M)0.5)+0; /生成单极性序列 else /双极性条件 a=sign(rand(1,M)0.5)-0.5); /生成双极性序列 end tmp=zeros(L,M); /一个码元的归零部分取零 L1=L*L0; /占空比，求出一个码元不归零部分的取样点数 tmp(1:L1,:)=ones(L1,1)*a; /将一个码元不归零部分的取样点值置为1 s=tmp(:); S=t2f(s,fs); /傅里叶变化 P=abs(S).2/T; /样本的功率谱密度 EP=EP*(1-1/loop)+P/loop; /随机过程的功率谱是各个样本的功率谱的数学期望end xset(window,1)plot(t,s)set(gca(),grid,1,1) title(时域图)xlabel(t)ylabel(S(t)mtlb_axis(-3,3,-1.5,1.5);xset(window,2)plot(f,abs(EP+%eps)set(gca(),grid,1,1) title(功率谱图形)xlabel(f)ylabel(功率) mtlb_axis(-35,35,-5,max(EP+%eps); xset(window,3) plot(f,10*log10(EP+%eps)set(gca(),grid,1,1) title(功率谱图形(dB)xlabel(f)ylabel(功率) 实验四：clear allexec t2f.sci ;exec f2t.sci ;/exec eyes.sci;function eyes(result,L,Na) N=length(result); tt=0:1:Na*L; set(gca(),auto_clear,off) for jj=1:Na*L:N-Na*L plot(tt,result(jj:jj+Na*L); endset(gca(),auto_clear,on)endfunctionk=input(请输入采样点数2k:（输入k值）);aa=0.5;N=2k;L=8;M=N/L;Rs=2; /码率2MbpsTs=1/Rs; /码元间隔dt=Ts/L; /时域采样间隔fs=1/dt;df=1/(N*dt); /频域采样间隔T=N*dt; /截断时间Bs=N*df/2; /系统带宽t=-T/2+dt/2:dt:T/2; /时域横坐标f=-Bs+df/2:df:Bs; /频域横坐标alpha=input(所需要的滚降系数是（01范围）：); / Hcos=zeros(1,N); i1=find(abs(f)(1-alpha)/(2*Ts)&abs(f)=(1+alpha)/(2*Ts); Hcos(i1)=Ts/2*(1+cos(%pi*Ts/alpha*(abs(f(i1)-(1-alpha)/(2*Ts); i2=find(abs(f)=(1-alpha)/(2*Ts); Hcos(i2)=Ts; / GT=sqrt(Hcos); GR=GT; / for loop1=1:20; Eb_N0(loop1)=(loop1-1); / eb_n0(loop1)=10(Eb_N0(loop1)/10); EP=1; n0=EP/eb_n0(loop1); / sita=n0*Bs; / n_err=0; / EP1=1;/zeros(1,N); EP2=1;/zeros(1,N); EP3=1;/zeros(1,N); for loop2=1:20; b=sign(rand(1,M,normal); s=zeros(1,N); / s(L/2:L:N)=b/dt; SS=t2f(s,fs); P=abs(SS).2/T; /序列功率谱密度 EP1=EP1*(1-1/loop2)+P/loop2; / S2=SS.*GT; / s2=real(f2t(S2,fs); /PAM P2=abs(S2).2/T; EP=EP*(1-1/loop2)+P2/loop2; / nr=sqrt(sita)*rand(1,N,normal); / sr=s2+nr; / SR=t2f(sr,fs); S=SR.*GR; PS=abs(S).2/T; / EP3=EP3*(1-1/loop2)+PS/loop2; / /a=real(f2t(S,fs); / sout=real(f2t(S,fs); y=sout(L/2:L:N); / bb=sign(y); / s3=zeros(1,N); / s3(L/2:L:N)=bb/dt; BB=t2f(s3,fs); PB=abs(BB).2/T; EP2=EP2*(1-1/loop2)+PB/loop2; / n_err=n_err+length(find(bb=b);/ end Pe(loop1)=n_err/(M*loop2); xset(window,10) plot(Eb_N0,log10(Pe+%eps),g); /PeEb/N0 xlabel(Eb/N0);ylabel(Pe);title(PeEb/N0曲线); eb_n0=10.(Eb_N0/10); set(gca(),auto_clear,off) plot(Eb_N0,log10(0.5*erfc(sqrt(eb_n0); mtlb_axis(0,15,-3.5,0); xlabel(Eb/N0) ylabel(Pe) legend(实际的,理论的); set(gca(),auto_clear,on) end xset(window,1) plot(f,EP1) title(序列功率谱) xlabel(f(kHz) ylabel(功率谱(W/kHz) mtlb_axis(-2,2,0,max(EP1) xgrid xset(window,2) plot(t,s) title(发送序列) xlabel(t(ms) ylabel(s(t)(V) mtlb_axis(0,5,-2,2) xgrid xset(window,3) plot(f,EP) title(发送PAM信号功率谱) xlabel(f(kHz