MATLAB通信系统仿真实验报告1

MATLAB通信系统仿真实验报告

实验一、MATLAB的基本使用与数学运算

目的：学习MATLAB的基本操作，实现简单的数学运算程序。

内容：

1-1要求在闭区间[0，2π]上产生具有10个等间距采样点的一维数组。试用两种不同的指令实现。

运行代码：x=[0:2*pi/9:2*pi]

运行结果：

1-2用M文件建立大矩阵x

x=[0.10.20.30.40.50.60.70.80.9

1.11.21.31.41.51.61.71.81.9

2.12.22.32.42.52.62.72.82.9

3.13.23.33.43.53.63.73.83.9]

代码：x=[0.10.20.30.40.50.60.70.80.9

1.11.21.31.41.51.61.71.81.9

2.12.22.32.42.52.62.72.82.9

3.13.23.33.43.53.63.73.83.9]

m_mat

运行结果：

1-3已知A=[5，6;7，8]，B=[9，10;11,12]，试用MATLAB分别计算

A+B,A*B,A.*B,A^3，A.^3，A/B,A\B.

代码：A=[56;78]B=[910;1112]x1=A+B X2=A-B X3=A*B X4=A.*B X5=A^3 X6=A.^3X7=A/B X8=A\B

运行结果：

1-4任意建立矩阵A，然后找出在[10,20]区间的元素位置。

程序代码及运行结果：

代码：A=[1252221417;111024030;552315865]c=A>=10&A<=20运行结果：

1-5总结：实验过程中，因为对软件太过生疏遇到了些许困难，不过最后通过查书与同学交流都解决了。例如第二题中，将文件保存在了D盘，而导致频频出错，最后发现必须保存在MATLAB文件之下才可以。第四题中，逻辑语言运用到了ij，也出现问题，虽然自己纠正了问题，却也不明白错在哪了，在老师的讲解下知道位置定位上不能用ij而应该用具体的整数。总之第一节实验收获颇多。

实验二、MATLAB程序的编写

目的：掌握顺序结构、选择结构、循环结构程序设计方法。学会编写函数。

内容：

2-1编写程序，建立向量N=[1,2,3,4,5]，然后利用向量N产生下列向量；（1）2,4,6,8,10

(2)1/2，1，3/2，2，5/2

(3)1，1/2，1/3，1/4，1/5

(4)1，1/4，1/9，1/16，1/25

代码：N=[1,2,3,4,5]X1=N*2X2=N/2X3=1./N X4=X3*X3

运行结果：

2-2从键盘输入一个三位整数，将他反向输出，如输入为639，输出为936.输入一个百分制成绩，要求输出成绩等级A,B,C,D,E。其中90~100分为A，80~89分为B，70~79分为C，60~69分为D，60分以下为E。

要求：

（1）分别用if语句

代码：clear

m=input('请输入一个三位数：')

m1=fix(m/100);

m2=rem(fix(m/10),10);

m3=rem(m,10);

n=m1+m2*10+m3*100;

disp(n);

（2）clear;

Mark=input('请输入成绩：');

Rank=cell(1,5);

S=struct('Marks',Mark,'Rank',Rank); for i=1:10;

a{i}=89+i;

b{i}=79+i;

c{i}=69+i;

d{i}=59+i;

e{i}=0+i;

q{i}=9+i;

g{i}=19+i;

h{i}=29+i;

m{i}=39+i;

n{i}=49+i;

end;

for i=1:5;

switch S(i).Marks

case100

S(i).Rank='A';

case a

S(i).Rank='A';

case b

S(i).Rank='B';

case c

S(i).Rank='C';

case d

S(i).Rank='D';

case e

S(i).Rank='E';

case q

S(i).Rank='E';

case g

S(i).Rank='E';

case h

S(i).Rank='E';

case m

S(i).Rank='E';

case n

S(i).Rank='E';

otherwise

S(i).Rank='成绩输入错误';

end

end

disp([num2str(S(i).Marks),blanks(3),S(i).Rank]);disp('');

运行结果：

2—3输入20个两位随机数，求其中的最大数最小数。要求分别用循环结构和调用MATLAB的max函数、min函数实现。

(1)a=fix(rand(1,20)*100)ma=max(a)mi=min(a)

运行结果：

（2）a=fix(rand(1,20)*100);

for i=1:20;

max=a(1);

min=a(1);

if max

max=a(i);

end

if min>a(i);

min=a(i);

end

end

max

min

运行结果：

2-6写出下列程序输出结果

(1)s=0;

a=[12,13,14;15,16,17;18,19,20;21,22,23]; for k=a

for j=1:4

if rem(k(j),2)~=0

s=s+k(j);

end

end

end

s

运行结果：

(2)global x

x=1:2:5;

y=2:2:6;

sub(y);

x

y

(3)function fun=sub(z)

global x

z=3*x;

x=x+z;

运行结果：

总结：第二次实验，对软件的使用比较熟练了，但还是遇到了些许问题。在运算符号的使用中，应当注意“.*”的使用，在最初因为不太会运用遇到了些困难，后来通过同学讨论和翻阅课本找到了答案。2—2中的第二种方法是按照课本例题改编的，有些啰嗦，不多至少是结果正确。还有2—6中刚开始没能正常输出，在老师的指导下知道（2）（3）是一起使用，算是运用到了函数调用。好在最后所有题目都得到了满意的结果。

实验三、MATLAB图形处理

目的：能够根据数据绘制各种形状的二、三维图形。

3-1绘制曲线y=x^3+x+1,x的取值范围为[-5,5]

代码：

x=-5:0.01:5

y=x.^3+x+1

plot(x,y)

运行结果：

3-4有一组测量数据满足y=exp(-a*t),t的变化范围为0~10，用不同的线性和标记点画出a=0.1,a=0.2和a=0.5三种情况下的曲线。

代码：

t=0:0.1:10;

y1=exp(-0.1*t);

y2=exp(-0.2*t);

y3=exp(-0.5*t);

title('t from0to10');

plot(t,y1,t,y2,t,y3);

xlabel('Variable t');

ylabel('Variable y');

text(0.8,1.5,'曲线y1=exp^{-0.1t}');

text(2.5,1.1,'曲线y1=exp^{-0.2t}');

text(0.8,1.5,'曲线y1=exp^{-0.5t}');

legend('y1','y2','y3')

运行结果：

3-7绘制饼图，x=[6649715638],并将第五个切块分离出来。

代码：

x=[6649715638];

subplot(1,2,1);

pie(x);

subplot(1,2,2);

pie(x,[0,0,0,0,1]);

运行结果：

总结：这次实验，比较有成就感，并没有遇到什么太复杂的困难，但是软件操作上出现了写麻烦，一不小心将软件页面的各个功能窗口关上了，颇费周折终于找到了那些功能窗口，但是整个页面都有些混乱。好在还是将题目做了出来，图出现的时候感觉特别有成就感。真的说明一件事情，英语学不好很麻烦啊。

实验四、MATLAB仿真模拟调制

目的：能用MATLAB仿真调幅信号和调角信号。

5-1用在区间[0，2]内的信号

m(t)=t0<=t<=1;m(t)=-t+21<=t<=2;

以DSB-AM方式调制一个载波频率为25HZ、幅度为1的载波产生已调信号u（t）。写一个Matlab的M文件，并用该文件作下面的题：

（1）画出已调信号；

（2）求已调信号的功率；

（3）求已调信号的振频谱，并与消息信号m(t)的频谱作比较。

程序代码：

dt=0.01;%时间采样间隔

fc=25;

T=1;

N=floor(T/dt);

t1=[0:N]*dt;

t2=t1+1;

%t=[t1t2];

mt1=t1;%信源

mt2=-t2+2;

%DSB-AM modulation

dsb1=mt1.*cos(2*pi*fc*t1);

dsb2=mt2.*cos(2*pi*fc*t2);

subplot(2,2,1);

plot(t1,dsb1);hold on;plot(t2,dsb2);

pwr1=mt1.^2;

pwr2=mt2.^2;

subplot(2,2,2);

plot(t1,pwr1);hold on;plot(t2,pwr2);

[mtf1,mtfft1]=FFT_SHIFT(t1,mt1);

[mtf2,mtfft2]=FFT_SHIFT(t2,mt2);

subplot(2,2,3);

plot(mtf1,abs(mtfft1));hold on;plot(mtf2,abs(mtfft2));

运行结果：

5-2设AM调整时，输入信号为没（t）=0.2sin1000pi*t+0.5cos1000exp2 *pi*t,A=1,载波中心频率fc=10khz

(1)用MATLAB画出AM信号的波形及其频谱

程序代码：

1、function[f,sf]=FFT_SHIFT(t,st)

df=t(2)-t(1);

T=t(end);

df=1/T;

N=length(t);

f=[-N/2:N/2-1]*df;

sf=fft(st);

sf=fftshift(sf);

2、dt=0.00001;%时间采样间隔

fm1=500;fm2=500*1.414;%信源频率

fc=10000;%载波中心频率

T=0.01;

N=floor(T/dt);

t=[0:N-1]*dt;

mt=0.2*sin(2*pi*fm1*t)+0.5*cos(2*pi*fm2*t);%信源

%AM modulation

A=1;

am=(A+mt).*cos(2*pi*fc*t);

[f,AMf]=FFT_SHIFT(t,am);

subplot(311);

plot(t,mt);

subplot(312);

plot(t,am);

subplot(313);

plot(f,AMf);

运行结果：

5-3设FM调制时，调频器的输入信号为一个周期性的锯齿波，锯齿波的一个周期为信号g(t)=t0<=t<1,g(t)=0其他，FM的中心频率fc=100hz,Kfm=10hz,试做（1）画出调频后的信号波形及其振幅谱

（2）若接收端采用鉴频器进行解调，且AWGN信道的功率密度谱为N0/2，试画出当解调器输入信噪比0dB,10Db,20dB时的解调输出信号，并与原信号进行比较。

程序代码：

1、function[f,sf]=FFT_SHIFT(t,st)

df=t(2)-t(1);

T=t(end);

df=1/T;

N=length(t);

f=[-N/2:N/2-1]*df;

sf=fft(st);

2.dt=0.001;%时间采样间隔

fc=100;

T=1;

N=floor(T/dt);

t=[0:N]*dt;

kf=10;

mt=t;

mti=t.^2/2;

fmt=cos(2*pi*fc*t+2*pi*kf*mti);

figure(1);

subplot(2,1,1);

plot(t,fmt);hold on;plot(t,mt,'r');

[f,ft]=FFT_SHIFT(t,fmt);

subplot(2,1,2);

plot(f,abs(ft));

s=1/2;%调制信号功率是A^2/2;

sn=100;db=s/(10^(sn/10));%求白噪声的方差。

noise1=sqrt(db)*randn(size(t));%产生高斯白噪声。fmt1=fmt+noise1;

figure(2);

subplot(2,1,1);

plot(t,fmt1);

N=length(fmt);

dfmt=zeros(1,N);

for k=1:N-1%已调信号微分

dfmt(k)=(fmt1(k+1)-fmt1(k))/dt;

end

envlp=abs(hilbert(dfmt));%求瞬时幅度，即包络subplot(2,1,2);

plot(t,envlp);

运行结果：

总结：这一次实验遇到的问题是花费了两节课解决的。第一个题目很快的完成后，第二个跟第三个题目都遇到了问题。第二题不出图，第三题图不正确，反复检查修改了好多遍程序还是不对，一直提示[f,ft]=FFT_SHIFT(t,fmt)程序这句话有问题，就在打算放弃的时候，终于在老师的提醒下知道，是因为没有调用函数，讲课本后面附录中的函数体加上之后，图就顺利的出来了。说明对程序的理解还不够到位，需要继续学习和努力。

实验五、MATLAB仿真模拟信号的数字传输

目的：能用MATLAB仿真函数的抽样、量化过程，掌握信号编码方法。6-1设低通信号s(t)=sin2*pi*t+0.5cos4*pi*t

(1)画出该低通信号的波形；

(2)画出抽样速率为fs=4Hz的抽样序列；

(3)从抽样序列恢复出原始信号；

(4)当抽样速率fs=2hz时，画出恢复出的抽样信号。

运行代码：1.clear all;close all;

dt=0.01;

t=0:dt:10;

xt=sin(2*pi*t)+0.5*cos(4*pi*t);

[f,xf]=FFT_SHIFT(t,xt);

fs=4;

sdt=1/fs;

t1=0:sdt:10;

st=sin(2*pi*t1)+0.5*cos(4*pi*t1);

[f1,sf]=FFT_SHIFT(t1,st);

t2=-50:dt:50;

gt=sinc(fs*t2);

stt=INSERT0(st,sdt/dt);

xt_t=conv(stt,gt);

figure(1)

subplot(3,1,1);

plot(t,xt);title('原始信号');

subplot(3,1,2);

stem(t1,st);title('抽样信号');

subplot(3,1,3);

t3=-50:dt:60+sdt-dt;

plot(t3,xt_t);title('抽样信号恢复');

axis([010-11])

2.function[out]=INSERT0(d,M)

N=length(d);

out=zeros(1,M*N);

for i=0:N-1

out(i*M+1)=d(i+1);

end;

3.function[f,sf]=FFT_SHIFT(t,st)

df=t(2)-t(1);

T=t(end);

df=1/T;

N=length(t);

f=[-N/2:N/2-1]*df;

sf=fft(st);

运行结果：

6-2用一个均匀量化器对零均值、单位方差的高斯源进行量化，这个量化器在区间[-10，10]内均匀量化。假定量化电平设在各量化区域的中间点，求出并画出量化电平数为N=3.4.5.6.7.8.9.10时，量化产生的均方失真作为量化电平数N 的函数。

t=[0:0.01:1];

s=normrnd(0,1,size(t));

err=zeros(1,8);

for i=1:8

[err(i),x_qtz]=myquantizer(s,i+2,-10.10);

end

N=[3,4,5,6,7,8,9,10];

Plot(n,err);

function[err,x_qtz]=myquantizer(x,n,l,h)

xmax=max(l,h);

x_qtz=x/xmax;delta=2/n;

q=delat*[0:n-1]-(n-1)/2*delat;

for i=1:n

index=find((q(i)-delat/2<=x_qtz)&(x_qtz<=q(i)+delta/2));

x_qtz(index)=q(i)*ones(1,length(index));end

x_qtz=x_qtz*xmax;

err=sum((x-x_qtz).^2)/length(x);

function[sqnr,x_qtz,code]=UniPcm(x,n)

xmax=max(abs(x));

x_qtz=x/xmax;

b_qtz=x_qtz;

delta=2/n;

q=delat*[0:n-1]-(n-1)/2*delta;

for i=1:n;

index=delat=find((q(i)-delta/2<=x_qtz)&x_qtz<=q(i)+delta/2));

x_qtz(index)=(q(i)*ones(1,length(index));

b_qtz(find(x_qtz==q(i)))=(i-1)*ones(1,length(find(x_qtz==q(i)))); end

x_qtz=x_qtz*xmax;

nu=ceil(log2(n));

code=zeros(length(x),nu);

for i=1:length(x)

for j=nu:-1:0

if(fix(b_qtz(i)/(2^j))==1)

code(i,nu-j)=1;

b_qtz(i)=b_qtz(i)-2^j;

end

end

end

6-6试编写A律13折线近似法的编码与解码程序。

function code=myAcode(x)

iw=[01632641282565121024]

实验六、MATLAB仿真数字信号的基带传输

目的：能够绘制常用码型，码型功率谱和眼图。

7-1设二进制符号序列为1011010010，以矩形脉冲为例，利用Matlab分别画出相应的单极性不归零、双极性不归零、单极性归零和双极性归零波形。

程序代码：

1、单极性不归零

function y=snrz(x);

t0=200;

t=0:1/t0:length(x);

for i=1:length(x)

if x(i)==1

for j=1:t0

y=((i-1)*t0+j)=1;

end

else

for j=1:t0

y=((i-1)*t0+j)=0;

end

end

end

y=[y,x(i)];

plot(t,y);

title(‘1011010010’);

axis([0,i,-0.1,1.1]);

2、双极性不归零：将程序1中的y=((i-1)*t0+j)=0改为y=((i-1)*t0+j)=-1

3、单极性归零

function y=srz(x);

t0=200;

t=0:1/t0:length(x);

for i=1:length(x)

if x(i)==1

for j=1:t0/2

y=((2*i-2)*t0/2+j)=1;

y=((2*i-1)*t0/2+j)=0;

end

else

for j=1:t0

y=((i-1)*t0+j)=0;

end

end

end

y=[y,x(i)];

plot(t,y);

title(‘1011010010’);

grid on;

axis([0,i,-0.1,1.1]);

4、双极性归零：将程序3中的for j=1:t0y=((i-1)*t0+j)=0改为

for j=1:t0/2

y=((2*i-2)*t0/2+j)=-1;

y=((2*i-1)*t0/2+j)=0;

运行结果：

7-4假设随机二进制序列“10110001”，“1”码对应的基带波形为升余弦波形，持续的时间为Ts;”0”码对应的基带波形与“1”码相反。画出基带波形以及眼图。

程序代码：

1、Ts=1;N=15;

eye_num=6;

a=1;N_data=1000;

dt=Ts/N,

t=-3*Ts:dt:3*Ts;

d=sign(randn(1,N_data));

dd=INSERT0(d,N);

ht=sinc(t/Ts).*(cos(a*pi*t/Ts))./(1-4*a^)*t.^2/Ts^2+eps);

st=conv(dd,ht);

tt=-3*Ts:dt:(N_data+3)*N*dt-dt;

subplot(211);plot(tt,st);axis([020-1.21.2]);

xlabel(‘t/Ts’);ylable(‘基带信号’);subplot(212);