通信原理编程实验完整代码

通信原理编程实验完整代码

告各位学弟学妹书：

此实验对于大多数人来说甚是烦人，但是我还是要劝你：

尽力去学吧，学的时候觉得没什么用，那是因为不了解其博

大精深，如果你一不小心顶住压力认真学了(譬如me),日

后涉及相关领域时，你就能够触类旁通，瞬间迸发出灵感，

撞击出智慧的火花，这种感觉就好像搂着一群漂亮的姑娘行

走在春风中一样，原来那些朦胧的理论都变成了清晰的线

条，串在一起：你会有一种类似于愚公般的给子子孙孙签订

卖身契投身于学术的激进想法，不过这种极左的思想毕竟要

不得一一除非你想女朋友和你分手，于激励学习便是极好

的。

实验一：

〃基带信号m(t)=sin(2000*pi*t)+2cos(1000*pi*t),fc=20khz,^<

AMZDSB-SC,SSB

clearall

exect2f.sci;

execf2t.sci;

fs=800;〃采样速率

T=200;〃截短时间

N=T*fs;〃采样点数

dt=l/fs;〃时域采样间隔

t=[-T/2:dt:T/2-dt];〃时域采样点

df=VT;〃频域采样间隔

f，・fs/2:df:fs/2-df];〃频域采样点数

fml=l;〃待观测正弦波频率，单位KHz,下同

fm2=0.5;〃待观测余弦波频率

fc=20;〃载波频率

〃以上为初始化参数设置

ml=sin((2*%pi)*fml*t);〃待观测正弦波部分

Ml=12f(ml,fs);〃傅里叶变换

MH1二・%i*sign(f).*Ml;〃希尔伯特变换

mhl=real(f2t(MHl,fs));〃希尔伯特反变换

m2=2*cos((2*%pi)*fm2*t);//待观测余弦波部分

M2=t2f(m2,fs);〃傅里叶变换

MH2=-%i*sign⑴.*M2;〃希尔伯特变换

mh2=real(f2t(MH2,fs));〃希尔伯特反变换

sl=(l+(ml+m2)/abs(max(ml+m2))).*cos((2*%pi)*fc*t);//AM

信号时域表达式

Sl=t2f(sl,fs);//AM信号频域表达式

s2=(ml+m2).*cos((2*%pi)*fc*t);//DSB-SC信号时域表达式

S2=t2f(s2,fs);//DSB-SC信号频域表达式

s3=(ml+m2).*cos((2*%pi)*fc*t)-(mhl+mh2).*sin((2*%pi)*fc*t);

〃SSB信号时域表达式，以上边带为例

S3=t2f(s3,fs);//SSB信号上边带频域表达式

〃以上是仿真计算部分

〃以下为绘图部分

//AM信号

xset('window',l)

plot(fzabs(Sl))

title('AM信号频谱')

xlabelCf)

ylabelCSff),)

mtlb_axis([-25,25,0,max(abs(Sl))]);

xsetCwindow',2)

plot(t,sl)

title('AM信号波形')

xlabel('t')

ylabel('s(t)')

；

mtlb_axis([-3,3,-3;3])

//DSB-SC信号

xset('window',3)

plot(fzabs(S2))

title('DSB-SC信号频谱')

xlabelCf')

ylabel('S(f)')

mtlbaxis([-25/25/0,max(abs(S2))]);

xset('window',4)

plot(t,s2)

title(DSB-SC信号波形)

xlabel('t')

ylabel('s(t)')

；

mtlb_axis([-l,4,-3;3])

〃SSB信号(以上边带为例)

xsetCwindow',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^),fc=

40KHz,Kf=5kHz/V,仿真产生FM,观察频谱和波形

〃初始化参数设置

clearall;

exect2f.sci;

execf2t.sci;

xdel(winsid());〃关闭所有图形窗口

fs=4000;〃采样速率

T=16;〃截短时间

N=T*fs;〃采样点数

dt=l/fs;〃时域采样间隔

t=[-T/2:dt:T2dt];〃时域采样点

df=l/T;〃频域采样间隔

f二1fs/2:df:fs/2-df];〃频域采样点数

fml=l;

fm2=0.5;

fm3=0.25;

fc=40;

Kf=5;〃频偏常数

〃仿真计算

ml=sin((2*%pi)*fml*t);

phil=(2*%pi)*Kf*cumsum(ml)*dt;//cumsum函数为求累积

和，相当于积分

m2=2*cos((2*%pi)*fm2*t);

phi2=(2*%pi)*Kf*cumsum(m2)*dt;

m3=4*sin((2*%pi)*fm3*t+%pi^);

phi3=(2*%pi)*Kf*cumsum(m3)*dt;

s=cos((2*%pi)*fc*t+phil+phi2+phi3);//FM信号时域表达式

S=t2f(s,fs);//FM信号频域表达式

m=sin(2*%pi*fml*t)+2*cos(2*%pi*fm2*t)+4*sin(2*%pi*fm3*t

+%pi3)；〃调制信号表达式

M=t2f(m,fs);

〃绘图

//FM已调信号

xset('window',l)

plot(f,abs(S))

titleCFM信号频谱")

xlabelCf')

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))

titleCm(t)信号频谱)

xlabel('f')

ylabel('S(f)')

mtlb_axis([-5,5,0/max(abs(M))]);

xset('window'M)

plot(t,m)

title(，m(t)信号波形)

xlabel('t')

ylabel('m(t)')

mtlbaxis([-8,8,-6,8]);

实验三：

〃通过仿真测量占空比为50%、75%以及100%的单双极性归

零码波形及其功率

〃谱，分析不同占空比对仿真结果的影响

clearall;

xdel(winsid());〃关闭所有图形窗口

exect2f.sci;

execf2t.sci;

L=32;〃每个码元间隔内的采样点数

N=2A13;〃总采样点数

M=N/L;〃总码元数

Rb=2;〃码元速率

Ts=VRb;〃比特间隔

fs=L/Ts;〃采样速率

T=N/fs;〃截短时间

Bs二fs〃；〃系统带宽

t=-T/2+[0:N-l]/fs;〃时域采样点

f=-Bs+[O:N-l]/T;〃频域采样点

L0二inputs请输入占空比(0~1)：')

EP=zeros(l,N);

ch=inpu"请选择要观察的码型：1-单极性；2.双极性)

forloop=l:1000//1000次样本函数取平均

ifch==l〃单极性条件

a=(rand(l,M)>0.5)+0;〃生成单极性序列

else〃双极性条件

3力811((1^或1,1\/1)>0.5)-0.5);〃生成双极性序歹1|

end

tmp=zeros(L,M);〃一个码元的归零部分取零

L1=L*LO;〃占空比，求出一个码元不归零部分的取样点

数

；〃将一个码元不归零部分的

tmp([l:Ll],:)=ones(Ll,l)*a

取样点值置为1

s=tmp(:)';

S工2f(s,fs);〃傅里叶变化

P=abs⑸.八如;〃样本的功率谱密度

EP=EP*(l-yioop)+P/loop;〃随机过程的功率谱是各个样

本的功率谱的数学期望

end

xset("window",l)

plotts)

set(gca()Jgrid"l,l])

title('时域图')

xlabel('t')

ylabegt),)

mtlb_axis([-3,3,-1.5,1.5]);

xset(,,window",2)

plot(f,abs(EP+%eps))

setfgcaO/'gridMl.l])

title。功率谱图形)

xlabel('f')

ylabel。功率)

mtlb_axis([・35,35u5,max(EP+%eps)]);

xset("window",3)

plot(f,10*logl0(EP+%eps))

set(gca()Jgrid"Ll])

title(,功率谱图形(dB),)

xlabel('f')

ylabel。功率)

实验四：

clearall

exect2f.sci;

execf2t.sci;

//execeyes.sci;

functioneyes(result,L,Na)

N=length(result);

tt=O:l:Na*L;

set(gca()/,auto_clear，7,off")

forjj=l:Na*L:N-Na*L

plot(tt,result(jj:jj+Na*L));

end

,

set(gca(),"auto_clear",'on")

endfunction

k=input(，请输入采样点数2八k:(输入k值));

aa=0.5;

N=2Ak;

L=8;

M=N/L;

Rs=2;〃码率2Mbps

Ts=VRs；〃码元间隔

dt=Ts/L;〃时域采样间隔

fs=]/dt;

df=l/(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('所需要的滚降系数是(0~1范围)：1);//

Hcos=zeros(l,N);

il=find(abs(f)>(l-alpha)/(2*Ts)&abs(f)<=(l+alpha)/(2*Ts));

Hcos(il)=Ts/2*(l+cos(%pi*Ts/alpha*(abs(f(il))-(l-alpha)/(2*Ts))

))；

i2=find(abs(f)<=(l-alpha)/(2*Ts));

Hcos(i2)=Ts;

//

GT=sqrt(Hcos);

GR=GT;//

forloopl=l:20;

Eb_NO(loopl)=(loopl-l);//

eb_n0(loopl)=10A(Eb_N0(loopl)/10);

EP=1;

nO=EP/eb_nO(loopl);//

sita=nO*Bs;//

n_err=0;//

EPl=l;//zeros(l,N);

EP2=l;//zeros(l,N);

EP3=l;//zeros(l,N);

forIoop2=l:20;

b=sign(rand(l,M；,normal"));

s=zeros(l,N);//

s(L/2:L:N)=b/dt;

SS=t2f(s,fs);

P二abs(SS).八NT;〃序列功率谱密度

EP1=EP1*(1-1/Ioop2)+P/Ioop2;//

S2=SS.*GT;//

s2=real(f2t(S2zfs));//PAM

R2=abs(S2).AyT;

EP=EP*(l-yioop2)+P^loop2;//

nr=sqrt(sita)*rand(l,N/,normal");//

sr=s2+nr;//

SR=t2f(sr,fs);

S=SR.*GR;

PS二abs⑸.八NT;//

EP3=EP3*(l-l/loop2)+PS/loop2;//

//a=real(f2t(S,fs));//

sout=real(f2t(S,fs));

y=sout(L/2:L:N);//

bb=sign(y);//

s3=zeros(l,N);//

s3(L/2:L:N)=bb/dt;

BB=t2f(s3,fs);

PB=abs(BB).A2/T;

EP2=EP2*(l-l/loop2)+PB/loop2;//

n_err=nerr+length(find(bb"=b));//

end

Pe(loopl)=n_err/(M*loop2);

xset("window",10)

plot(Eb_NOJoglO(Pe+%eps)/g');//Pe^Eb/NO

xlabelCEb/NO'hylabelCPe^jtitleC'Pe-Eb/NO曲线)

eb_n0=10.A(Eb_N0/10);

set(gca()/'auto_clear","off")

plot(Eb_N0,logl0(0.5*erfc(sqrt(eb_n0))));

mtlb_axis([0,15,-3.5,0]);

xlabel('Eb/NO')

ylabel('Pe')

legend]实际的理论的');

,,

set(gca(),"auto_clear"/on")

end

xsetC'window",!)

plot(f,EPl)

title(”序列功率谱”)

xlabel("f(kHz)")

ylabel(“功率谱(W/kHz)”)

mtlb_axis([-2,2/0,max(EPl)])

xgrid

xset("window",2)

plot(t,s)

title("发送序列“)

xlabel("t(ms)")

ylabel(ns(t)(V)")

mtlb_axis([0,5,-2,2])

xgrid

xset("window",3)

plot(f,EP)

title("发送PAM信号功率谱”)

xla