南京邮电大学matlab软件设计超详细.doc

南京邮电大学matlab软件设计超详细.doc

南京邮电大学通信学院

软件课程设计
实验报告

模块名称：___MATLAB软件设计
专业班级：__通信工程
姓 名：____ ____
学 号：___

实验日期： 6 月 17—28日 实验报告日期： 年 7 月 1 日
一、要求练习的实验部分
1． 在时间区间 [0,10]中，绘出曲线。
程序：
t=0:0.1:10;
y=1-exp((-0.5)*t).*cos(2*t);
plot(t,y,'r-');
shg
结果：

2. 写出生成如图E2-1所示波形的MATLAB脚本M文件。图中虚线为正弦波，要求它的负半波被置零，且在处被削顶。
程序：
t=linspace(0,3*pi,500);
y=sin(t);
a=sin(pi/3);
z=(y>=0).*y;
z=(y>=a).*a+(y<a).*z;
plot(t,y,':r');
hold on;
plot(t,z,'b-')
xlabel('t'),ylabel('z=f(t)'),title('逐段解析函数')
legend('y=sin(t)','z=f(t)',4),hold off
结果：

3. 令，运行 [U,S,V]=svd(A);tol=S(1,1)*3*eps; 然后回答以下问题：
（1） sum(diag(S)>tol) - rank(A) 是多少 ？
（2） S(1,1) - norm(A) = 0 是多少 ?
（3） sqrt(sum(diag(S*S))) - norm(A,'fro') 的结果是什么 ?
（4） S(1,1)/S(3,3) - cond(A) 的结果是什么 ?
（5） S(1,1)*S(2,2)*S(3,3) - det(A) < tol 的运行结果是什么？
（6） V(:,1)'*null(A) 得到什么结果 ?
（7） abs(A*null(A)) < tol 得到什么结果 ?
（8） U(:,1:2) = = orth(A) 的运行结果是什么 ?

程序：

clear;
clc;
disp('设 A=[1,2,3;4,5,6;7,8,9]，得')
A=[1,2,3;4,5,6;7,8,9] %一个3*3矩阵
disp('设 [U,S,V]=svd(A)，得')
[U,S,V]=svd(A) %奇异值分解(SVD)USV:U m*m酉矩阵 V n*n酉矩阵 S 对角阵
disp('设 tol=S(1,1)*3*eps，得')
tol=S(1,1)*3*eps
disp('设 a=sum(diag(S)>tol) - rank(A)，求 a')
a=sum(diag(S)>tol) - rank(A)
disp('设 b=S(1,1) - norm(A)，求 b')
b=S(1,1) - norm(A)
disp('设 c=sqrt(sum(diag(S*S))) - norm(A, fro )，求 c')
c=sqrt(sum(diag(S*S))) - norm(A,'fro')
disp('设 d=S(1,1)/S(3,3) - cond(A)，求 d')
d=S(1,1)/S(3,3) - cond(A)
disp('设 e=S(1,1)*S(2,2)*S(3,3) - det(A) < tol ，求 e')
e=S(1,1)*S(2,2)*S(3,3) - det(A) < tol
disp('设 f=V(:,1) *null(A)，求 f')
f=V(:,1)'*null(A)
disp('设 g=abs(A*null(A)) < tol，求 g')
g=abs(A*null(A)) < tol
disp('设 h=U(:,1:2) == orth(A)，求 h')
h=U(:,1:2) == orth(A)
结果：
设 A=[1,2,3;4,5,6;7,8,9]，得

A =

1 2 3
4 5 6
7 8 9

设 [U,S,V]=svd(A)，得

U =

-0.2148 0.8872 0.4082
-0.5206 0.2496 -0.8165
-0.8263 -0.3879 0.4082

S =

16.8481 0 0
0 1.0684 0
0 0 0.0000

V =

-0.4797 -0.7767 -0.4082
-0.5724 -0.0757 0.8165
-0.6651 0.6253 -0.4082

设 tol=S(1,1)*3*eps，得

tol =

1.1223e-014

设 a=sum(diag(S)>tol) - rank(A)，求 a

a =

0

设 b=S(1,1) - norm(A)，求 b

b =

0

设 c=sqrt(sum(diag(S*S))) - norm(A, fro )，求 c

c =

3.5527e-015

设 d=S(1,1)/S(3,3) - cond(A)，求 d

d =

-8

设 e=S(1,1)*S(2,2)*S(3,3) - det(A) < tol ，求 e

e =

1

设 f=V(:,1) *null(A)，求 f

f =

0

设 g=abs(A*null(A)) < tol，求 g

g =

1
1
1

设 h=U(:,1:2) == orth(A)，求 h

h =

1 1
1 1
1 1
结果分析：
由上面的结果得出如下的题目答案：
（1）sum(diag(S)>tol) - rank(A) 的结果是0；
（2）S(1,1) - norm(A) = 0 的结果是0；
（3）sqrt(sum(diag(S*S))) - norm(A,'fro') 的结果是3.5527e-015；
（4）S(1,1)/S(3,3) - cond(A) 的结果是-8；
（5）S(1,1)*S(2,2)*S(3,3) - det(A) < tol 的运行结果是1；
（6）V(:,1)'*null(A)的结果是0；
（7）abs(A*null(A)) < tol的结果是1
1
1 ；
（8）U(:,1:2)==orth(A)的运行结果是1 1
1 1
1 1
4. 求积分，。
程序：

t=0:pi/100:2*pi;
f=abs(sin(cos(t)));
F=cumsum(f)*pi/100;
plot(t,F);
xlabel('x轴');
ylabel('y轴');
title('描绘y=abs(sin(cos(t)))在(0,x)间积分曲线[x∈(0,2π)]')
grid on;
结果：

5. 求方程的解。
程序：
clear;
clc;
disp('计算方程组 x^2+y^2=1 x*y=2 的根 x y')
[x,y] = solve('x^2+y^2 = 1','x*y = 2')
结果：
计算方程组 x^2+y^2=1 x*y=2 的根 x y

x =

[ -1/2*(1/2*5^(1/2)+1/2*i*3^(1/2))^3+1/4*5^(1/2)+1/4*i*3^(1/2)]
[ -1/2*(1/2*5^(1/2)-1/2*i*3^(1/2))^3+1/4*5^(1/2)-1/4*i*3^(1/2)]
[ -1/2*(-1/2*5^(1/2)+1/2*i*3^(1/2))^3-1/4*5^(1/2)+1/4*i*3^(1/2)]
[ -1/2*(-1/2*5^(1/2)-1/2*i*3^(1/2))^3-1/4*5^(1/2)-1/4*i*3^(1/2)]

y =

[ 1/2*5^(1/2)+1/2*i*3^(1/2)]
[ 1/2*5^(1/2)-1/2*i*3^(1/2)]
[ -1/2*5^(1/2)+1/2*i*3^(1/2)]
[ -1/2*5^(1/2)-1/2*i*3^(1/2)]
6. 在某激励条件下，二阶系统归一化响应可表示为，其中为阻尼系数，，。请用不同的颜色或线型，在同一张图上，绘制取值下系统在区间内的响应曲线，并要求用和对它们相应的两条曲线进行醒目的文字标志。
程序：
b=0.2:0.2:2.0;
t=0:0.1:18;
color=['r','g.','bo','cx','m+','y*','kv','r--','g:','b-.'];
for i=1:10
p(i)=sqrt(1-b(i).*b(i));
q(i)=atan(sqrt(1-b(i).*b(i))./b(i));
y=1-1./p(i).*exp(-b(i).*t).*sin(p(i).*t+q(i));
if i==1
plot(t,y,color(i))
text(t(30),y(30),' \leftarrow ξ=0.2','FontSize',18)
hold on
else if i==10
plot(t,y,color(i))
text(t(50),y(50),' \leftarrow ξ=2.0','FontSize',18)
hold on
else
plot(t,y,color(i))
hold on
end
end
end
结果：

7. 构建一个简单的全波整流模型，并用示波器分两路同时观察原信号和整流后的信号波形。要求：信源模块被重新命名为“输入正弦波”；信宿模块被重新命名为“示波器”；连接到信宿的信号线上分别标注“原信号”和“整流信号”；在模型窗中添加注释文本。
(1)、原理图：
各模块参量设置：
[示波器]
Parameters
Number of axes: 2
其它默认
[其余模块]
默认选项

结果

(2)、原理图：
各模块参数设定
[Gain]
Gain: -1
其它默认
[其余模块]
默认选项

结果

8. 利用 SIMULINK及其标准模块设计一个低通滤波器，从受噪声干扰的多频率混合信号中获取10Hz的信号。在此，而
各模块参量设置：
[Sine Wave]
Frequency (rad/sec): 10
其它默认
[Sine Wave1]
Frequency (rad/sec): 100
Phase (rad): pi/2
其它默认
[Random Number]
Variance: 0.2
Sample time: 100
其它默认
[Sum]
List of signs: |+++
其它默认
[Analog Filter Design]
Passband edge frequency (rad/s): 10
其它默认
[其余模块]
默认选项

结果：

9. 已知系统的状态方程为：，其中，请构建该系统的仿真模型，并用XY Graph模块观察相轨迹。
原理图：
参数配置;
[Sum]
List of signs: |--+
其它默认
[integrator]
Initial condition:0.2
[integrator1]
Initial condition:0.2

结果：

10. CD74HC00仿真（数字电路）

实现方法：
这里就是四个与非门的实现，利用Logical Operator模块做出四个与非门，还要再和Vcc和Gnd相与，再输出。然后用Subsystem封装成子系统即可；
原理图：
(1)、封装前的原理图

（2）、分装后的测试原理图

结果：

分析：前两路分别为产生“1100”和“0101”的脉冲发生器的波形，后四路波形则分别为四个与非门的输出；
能够看到，四路输出信号与输入信号均符合与非逻辑的关系，符合7400的真值表，设计正确！
11．2FSK信号的生成与接收滤波
要求：
1) 产生二进制[0，1]随机信号波，信号波特率为1000B
2) 对此信号进行2FSK调制，f1=1500Hz, f2=3050Hz
3) 加入高斯白噪声
4) 在接收端进行带通滤波
5) 抽样速率为20KHz.
6) 观察滤波前后信号的波形和频谱。
实现方法：
利用Bernoulli Binary Generator模块产生随机信号；
利用Sine Wave模块产生f1=1500Hz, f2=3050Hz 的载波，结合Switch模块进行2FSK调制；
利用AWGN Channel模块构建高斯白噪声信道，其参数设置为：信噪比：10dB，输入信号功率：1W；
利用Digital Filter Design设计带通滤波器，分别为1500Hz和3050Hz滤波
然后利用相干解调，在将两个载波相加，用符号函数判决出原信号。
原理图：（zoushuSimulink11.mdl、zoushuSimulink11b.mdl）
参数设计：
[Sine Wave]:
Frequency:1500 Hz
其它：默认
[Sine Wave 1]:
Frequency:3050 Hz
其它：默认
[Bernoulli Binary Generator]
Sample time:1/1000
其它：默认
[Switch]:
Threshoid:1/2
其它：默认
[AWGN Channel]
SNR:10 db
Power:1W
其它：默认
[Digital Filter Design ]
Bandpass;
Fs:20kHz
Fstop1:0.8KHz
Fpass1:1.2KHz:
Fstop2:2KHz
Fpass:2.2.5KHz:
其它：默认
[Digital Filter Design 1 ]
Bandpass;
Fs:20kHz
Fstop1:2.5KHz
Fpass2:2.8KHz:
Fstop2:3.5KHz
Fpass2:4KHz:
其它：默认
调制的原理图：

调制加解调的原理图：

结果：
频谱图：

调制、解调的波形图：

结果分析：原信号经过2FSK调制，经过噪声信道，再进行滤波、相干解调，最后又恢复出原信号，结果符合题目要求，系统设计正确！滤波器等相关设备的参数也设置恰当！实验成功！
12. 创立一个简单的离散多速率系统：单位阶跃信号经过具有不同速率的采样后分别用作两个离散传递函数的输入。这两个离散传递函数有相同的有理分式，但采样时间和时间偏置二元对分别设为[1 , 0.1]和[0.7 , 0]。要求：观察这两个离散传递函数的输出有什么不同；用不同的颜色标帜不同采样速率系统。（通信系统）

原理图：

结果

二、模拟数字电路仿真实验
(一) 数字逻辑电路基础
熟悉常见逻辑单元的特性，学会运用基本逻辑单元（与、或、非、异或、R-S触发器、D触发器、J-K触发器等），修改参数、增减输入/输出端。
(二) 组合逻辑电路仿真
1、 设计二/四线译码器
实现方法：
利用Logic Operator模块实现非门和与非门；
利用Pulse Generator模块产生0101和0011序列码，用来检验设计出的子系统。
原理图：

封装前原理图：

结果及分析

分析：当!EN端置0时，此时，随着第二、三两路信号输入的不同，作为输出端的Y0,Y1,Y2,Y3依次输出低电平，符合真值表要求，证明设计的子系统工作正常。

分析：当!EN端置1时，子系统没有被使能。此时，作为输出端的Y0,Y1,Y2,Y3没有随着输入A0，A1作出相应变化，而是一直输出高电平，证明设计的子系统工作正常。设计正确！

2、 设计四选一数据选择器
能仿真测试，并设计成子系统元件
实现方法：
利用的模块同上；
地址信号用Pulse Generator产生；
D0到D3也由Pulse Generator模块产生，为了区分各路信号，将其信号周期依次设置为0.0625s, 0.125s, 0.25s和0.5秒。

原理图

结果及分析：

前两路分别为输入信号，组合产生00，01，10，11的地址信号；第四到第七这四路信号为D0到D1，它们的频率各不相同；最后一路为输出信号；
能够看出，随着输入的地址信号的变化，输入一次输出该地址信号对应的信号。子系统设计成功！！
(三) 时序逻辑电路仿真
1、 设计四位二进制计数器（带置位和清零）
实现方法：
利用JK触发器和逻辑门构成之。利用比较系数的方法可得电路的驱动方程为

检验电路将P, T, LD, CR置为1，然后观察示波器输出。

原理图：
封装前的原理图：

封装后：

结果：

从图中能够看出，脉冲在上升沿有效，计数器按0000~1111方式进行计数。
当P、T、L（置数）、CR（清零）有一个为低电平时的输出为全零。P、T为0时，计数器不能工作；L为0置数，D3~D0全零，输出也全零。CR为零时清零，因此输出也全零。

三、数字信号处理仿真实验
（一）、利用Kaiser窗函数设计的FIR低通滤波器进行数字滤波
利用Kaiser窗函数，设计具有如下指标的FIR低通滤波器：fs=20KHz, fpass=4kHz, fstop= 5KHz, Apass=0.1dB, Astop=80dB。

归一化低通滤波器的主要性能指标有：
绝对指标：[0, wp]为滤波器的通带，δp为可容许的通带波动；[ws,π]为滤波器的阻带，δs为阻带波动。
相对指标：Ap为用dB表示的通带波动；As为用dB表示的阻带波动
相对指标与绝对指标的关系为：
Ap=-20lg((1-δp)/(1+δp))
As=-20lg(δs)
设计基本思路：首先选择一个符合要求的理想滤波器（这里是理想低通滤波器）；由于理想滤波器的冲击响应是非因果的且无限长，为了能用FIR滤波器实现，必须用适当的窗函数来截取，从而得到线性相位和因果的FIR滤波器。而滤波器的特性与窗函数的长度（也就决定了FIR滤波器的长度）密切相关。
一个截止频率为wc的理想低通滤波器的单位冲击响应为：
h(n)=sin[wc(n-a)]/[π(n-a)]
a为采样延迟，对应的频谱特性为：H(w)=exp(-jaw), 当|w|<wc时。
用一个窗函数w(n)去截取后，得到：
hd(n)=h(n).w(n)
窗函数满足在[0，M]内关于a对称，在其它区域为0。在频域，FIR滤波器的频率响应由理想滤波器的频响和窗函数的频响卷积得到。
对给定的阻带衰减，Kaiser窗提供了最大的主瓣宽度，从而提供最陡的过渡带。窗函数的表示式为：
w(n)=I0{βsquare(1-(1-2n/M)2)}
其中I0为零阶Bessel函数。Kaiser窗的优点在于它能够经过改变参数β和M来改变阻带衰减和过渡带宽。在MATLAB的SimLink中，有相关的模块供调用。
题目：
1 生成2KHz和8KHz的混合信号，使该信号经过上述滤波器，观察滤波前后的波形变化和频谱分布情况。（基本要求）
实现方法：首先用两个Sine Wave模块分别产生2kHz和8kHz的正弦信号，之后用add模块将二者混合。用Digital Filter Design设计低通滤波器，具体设置如图所示：
原理图：、

结果、

图为滤波之后和之前的信号频谱，能够观察到，滤波之后，处于通带之外的8kHz的正弦分量消失了。

上图是示波器的波形，其中第一路为2kHz和8kHz的混合信号，第二路为滤波之后的信号，可见低通滤波器符合要求！
2 在改信号中加入高斯白噪声，观察滤波前后的波形变化和频谱分布情况。（基本要求）
实现方法：在上题基础上增加AWGN信道，加入高斯白噪声。
原理图：、

结果：、

分析：图为滤波之后和之前的信号频谱，能够观察到，滤波之后，处于通带之外的8kHz的正弦分量消失了。

分析：上图是加噪声后示波器的波形，其中第一路为2kHz和8kHz和噪声的混合信号，第二路为滤波之后的信号（还有低频的噪声）！

3 、将上述滤波器改为中心频率为100KHz的带通滤波器，信号源为带外的90KHz和带内的100KHz的混合正弦信号，重做1，2问。
实现方法：
利用Digital Filter Design设计带通滤波器，具体设置为：Fs=1000kHz, Fstop1=99 kHz, Fpass1=99.5 kHz, Fpass2=100.5 kHz, Fstop2=101 kHz。

原理图：、

结果：
无噪声的频谱分析：

无噪声的波形分析：

加入噪声的频谱分析：

加入噪声的波形分析：

（二）、 DSB-SC信号的生成与解调
要求：
1) 用离散（DSP）的方法生成DSB信号
2) 载波频率为150KHz，音频为500Hz和 Hz的混合音。
3) 加入高斯白噪声
4) 带通滤波，分别用正交解调和检波的方法实现解调。
5) 抽样速率为600Hz。
6) 观察各点的信号波形和频谱。
A.正交解调法
实现方法
本题目要求用离散的方法生成DSB信号，因此我们用编写M文件方式而不是用Simulink来实现。生成正弦音频信号时，对时间t用向量表示，能够实现离散化。DSB调制与解调分别使用amod,ademod函数。加入高斯白噪声用awgn函数。带通滤波用cheb1ord和fftfilt函数实现。画频谱用fft函数实现。
程序代码（zoushuCode3002.m）：
t=0:1/(600*1000-1):0.01;
m1=sin(500*2*pi*t);%500Hz正弦信号
m2=sin(2*pi* *t);% Hz正弦信号
m=m1+m2;%混合信号
figure(1);
subplot(3,1,1);
plot(t,m);
title('原始信号');
z=amod (m,1.5*10^5,6*10^5);%调制
subplot(3,1,2),plot(t,z);
title('DSB-SC信号图');
y=awgn(z,5); %加入高斯白噪声
subplot(3,1,3),plot(t,y);
title('加入噪声后的DSB-SC图，SNR=5');
figure(2);
f=(0:length(fft(z))-1)*(600*1000-1)/length(fft(z));
subplot(3,1,1);plot(f,abs(fft(z)));
title('调制后的频谱图');
%带通滤波器参数
As=100; Ap=1;
fs=6*10^5; %采样频率
fp1=400; fp2=1.6*10^5; %通频带
fs1=300; fs2=1.7*10^5; %截止带
wp=[fp1 fp2]/(fs/2);
ws=[fs1 fs2]/(fs/2);
[n,wn]=cheb1ord(wp,ws,Ap,As);
b=fir1(n,wn,chebwin(n+1,As));
x=fftfilt(b,y);
z1=ademod(x,1.5*10^5,6*10^5);%滤波法解调
subplot(3,1,2);plot(t,z1);
title('解调后的波形')
f=(0:length(fft(z1))-1)*(600*1000-1)/length(fft(z1));
subplot(3,1,3);
plot(f,abs(fft(z1)));
axis([0 6000 0 3500]);
title('解调后的频谱图');

结果及分析：

加入高斯白噪声后对源信号干扰比较大，解调出来的波形失真较严重，但解调后的频谱图是正确的。
B.检波法
实现方法：
此方法与正交解调方法类似，只不过在解调时加入了直流信号，以便取出包络，取出包络后再减去直流信号得到源信号。取包络用取绝对值函数abs实现。
程序：（zoushuCode300202.m）
A=3;%包络解调加入直流信号
t=0:1/(600*1000-1):0.01;
m1=sin(500*2*pi*t);%500Hz正弦信号
m2=sin(2*pi* *t);% Hz正弦信号
m=m1+m2+A;%包络解调加入直流信号
figure(1);
subplot(3,1,1);%混合信号
plot(t,m);
title('和直流信号叠加在一起的原始信号');
z=amod(m,1.5*10^5,6*10^5,'amdsb-sc');%调制
subplot(3,1,2),plot(t,z);
title('DSB-SC信号图');
y=awgn(z,5); %加入高斯白噪声
subplot(3,1,3),plot(t,y);
title('加入噪声后的DSB-SC图，SNR=5');
figure(2);
f=(0:length(fft(z))-1)*(600*1000-1)/length(fft(z));
subplot(3,1,1);plot(f,abs(fft(z)));
title('调制后的频谱图');
%带通滤波器参数
As=100; Ap=1;
fs=6*10^5; %采样频率
fp1=400; fp2=1.6*10^5; %通频带
fs1=300; fs2=1.7*10^5; %截止带
wp=[fp1 fp2]/(fs/2);
ws=[fs1 fs2]/(fs/2);
[n,wn]=cheb1ord(wp,ws,Ap,As);
b=fir1(n,wn,chebwin(n+1,As));%fir滤波器
z1=pi/2*fftfilt(b,abs(y));%找出包络并进行缩放
subplot(3,1,2);plot(t,z1);
title('包络图')
z1=z1-A;%去除直流信号
f=(0:length(fft(z1))-1)*(600*1000-1)/length(fft(z1));
subplot(3,1,3);
plot(f,abs(fft(z1)));
axis([0 6000 0 3500]);
title('解调后的频谱图');
实验结果及分析

四、通信系统实验
1． 模拟调制---DSB调制与解调
要求：分别以
（1） 正弦信号（300Hz）
（2） 语音信号（自行录制’.wav’文件）
进行DSB调制和解调。信道采用高斯白噪声信道。
（1）正弦信号（300Hz）

结果：（图图4001.jpg）

（2）语音信号（自行录制’.wav’文件）
实现方法：
自己录制xy.wav文件放在E盘，用wavread读取语音文件。具体做法为在命令窗输入：y=wavread('F:\xy.wav',10001);t=[0:0.001:10]';运行，等到y为[10001*2],t为[10001*1]的数据。在simulink里调用From Workspace模块，将其Data属性设为：[t y]。用一个10KHz的正弦信号相乘去调制和解调。

结果分析：

结果分析：经过录制的声音，进过调制和解调，最后经过滤波，滤出原声音信号，可是由于录制的声音音量过小，最后滤出的波形没有那么精确，出现一混频。

五、学习心得：
为期两周的软件设计落下帷幕，这两种我从仅仅懂一些matlab皮毛知识到自己能够独立完成设计滤波器与通信系统与数字模块，之间也经历了很多挫折。其中最大的原因是在于自己对所要设计的原理不清楚，为此我重新捡起了数字电路、通信原理的书，从原理一点一滴做起，遇到实在不懂得，就去请教老师同学与百度，在于别人交流的过程中自己的设计思维也更明晰，少走了很多弯路。
例举一个遇到的小问题，在做最后一题通信仿真时，刚开始，在进行WAV格式语音信号处理的时候，程序始终现实如下错误提示：
    ??? Error using ==> wavread
    Error using ==> wavread
    Data compression format (IMA ADPCM) is not supported.
我在查阅了很多资料，在网上也查阅相关信息，花费了大量时间也没找出结果，最后发现在WAV格式的语音文件有两种格式，即PCM格式和IMA  ADPCM格式，而在MATLAB中用wavread函数进行语音处理时，并不能直接处理IMA ADPCM格式的语音信号，经过格式转换之后（选择PCM格式），我们运行出了正确的结果。实验过程中，我很清楚的感觉到初始语音信号和滤波输出后的语音信号在音色上有一定的差别，这说明了信号在处理、传输过程中有损耗。
软件设计的两周使我完成了从纯理论到仿真实践的第一步，这也是走向产品设计的很重要一步，使我对自己所学知识有了很好的巩固，总之，本次软件设计周很有意义。