数字中频GPS信号的MATLAB仿真1
杨勇，陈偲，王可东
北京航空航天大学宇航学院，北京 (100083)
E-mail：wangkd@https://www.doczj.com/doc/bd2793534.html,
摘要：文章以INS/GPS紧耦合为应用对象，在分析中频GPS信号结构的基础上，根据实际环境和载体运行状态，给出GPS信号延时、多普勒频移和钟差等参数，并应用中频信号解析表达式实现多颗卫星信号的合成。最后，基于MATLAB语言进行了仿真计算，仿真结果表明信号符合实际情况，同时经过软件接收机的捕获、跟踪和解调计算，验证了信号的正确性。
关键词：GPS；高动态；紧耦合；中频；信号模拟器
中图分类号：TP391
1．引言
随着固体弹道导弹射程的不断增加和打击精度的要求提高，纯惯性导航早已不能够满足要求。全球定位系统（GPS）和惯性系统（INS）相结合是复合制导的重要发展方向之一，而对于GPS/INS组合导航来说，为了缩短研制周期，便于新信号开发及测试，软件信号模拟器和接收机的研究成为重要的研究方向之一。GPS技术成长非常迅速，现在市场上的手持式GPS接收机已相当普遍，但是国内的自主知识产权的GPS技术产品的研发仍然比较薄弱，尤其是核心芯片的知识产权很少被国内所拥有。国内的“北斗”、“GALILEO”导航定位都处在发展之中，信号模拟器的研究被越来越多的被重视。
信号模拟器具有成本低、可重复性好、数据完整等优点，不仅能用于组合导航技术研究，也能为新信号的验证研究提供支持，还可以为硬件接收机的接收性能测试提供有效的信号环境模拟。
GPS信号模拟器是软件无线电研究的一个方面，为处于设计阶段的GPS接收机提供仿真环境。常见的GPS信号生成器产生的是射频信号，而目前接收机的设计重点侧重于基带信号处理，也就是本文中提到的数字中频GPS信号。数字中频GPS信号模拟器目前主要是仿真载体运动、模拟时钟偏差、卫星星钟误差、电离层误差、对流层误差、多路经效应、天线的方向、弹体振动、以及云层、雷雨等实际环境对GPS信号的影响，并对接收机前端的下变频、滤波、采样和自动增益控制进行仿真，直到生成GPS接收机信号处理所需的数字信号。
目前，国外已成功地开发出多种信号模拟器[1]；国内对高动态GPS信号的研究也比较深入，完成了星历的生成、动态信号的原理研究、误差建模和信号的仿真验证等研究工作 [2,3]。但国内研究的重点集中在数字中频GPS信号的生成，即针对环境误差、载体运动以及卫星星历的仿真，从而通过接收机捕获、跟踪、解算获得所需要的载体位置。国内大部分研究都成功仿真了GPS信号的功率谱，但是对于其信号真正用于导航解算的介绍不多，尤其是信号的实时性问题。本文主要针对紧耦合导航中对GPS信号生成的要求，分析GPS信号的生成与应用问题，并结合中频GPS信号的解析表达式，通过解算GPS信号的延迟，基于MATLAB语言，对大机动条件下的GPS信号进行仿真，获得紧耦合导航所需中频信号。
1本课题得到高等学校博士学科点专项科研基金（项目编号：20070006006）的资助。
2．GPS 中频信号的MATLAB 实现
2.1 GPS 中频信号
GPS 信号为了避免受到干扰，采用高频信号作为载波，对伪随机码和数据码进行调制，本文仅考虑使用L1载波（1575.42MHz ）的GPS 信号。由于L1频率很高，直接进行数字仿真在计算速度和精度上都无法保证，所以，这里直接对信号的中频解析表达式进行仿真，如式（1）所示。
}
IF tropo iono SV MP 1
tropo iono SV MP IF IF 1tropo iono SV MP 0()
()
cos[()()()]N
i p i i p d L p S t T t t t t C t T t t t t t t t T t t t t δδδδδδδδωωωωδδδδ?==???+????+?+?+???+??++∑（1）
其中：
r P 为信号功率，单位w ； ()i D 为数据码，频率50Hz ； ()i C 为伪随机码，频率1.023MHz ；
IF ω为中频信号角频率，单位rad/s ； d ω为多普勒频移角频率，单位rad/s ； 1L ω为L1载波的角频率，单位rad/s ； 0?为初始相位，单位rad ；
t 为接收机时刻，单位s ；
t ?是t 的整毫秒的小数部分，单位s （例如t =105.1121ms,t ?=0.1121ms ）
； p T 为从卫星发射到接收机接收的传输时间（不考虑误差），单位s ；
tropo t δ为对流层延时，单位s ； iono t δ为电离层延时，单位s ； SV t δ为星钟误差，单位s ； MP t δ为多路径效应误差，单位s ；
N 表示接收到的卫星数目。
为了获得所需的GPS 中频信号，同时简化问题，假设： （1）接收机钟差为零； （2）星钟误差SV 0t δ=； （3）多路径效应误差MP 0t δ=。 简化之后获得式（2）：
}
IF tropo iono 1
tropo iono IF IF 1tropo iono 0()
()
cos[()()()]N
i p i i p d L p S t T t t C t T t t t t t T t t δδδδωωωωδδ?==??????+?+???+?+∑
（2）
根据式（2）可得到最基本的载波信号的关系。在实际应用中，该式也能应用于其他简单模型。例如考虑了星钟误差之后，可利用式（1）进行仿真；如果考虑因水面折射、建筑物反射等原因引起的多路径效应时，可用式（3）表示。其中i α表示信号功率增益系数。
}{
IF tropo iono 1
tropo iono IF IF 1tropo iono 0tropo iono MP 1tropo iono MP IF IF ()
()
cos[()()()]()
()
cos[()()N
i p i i p d L p N
i p i i p d S t T t t C t T t t t t t T t t t T t t t C t T t t t t t t δδδδωωωωδδ?αδδδδδδωωω===??????+?+???+?++????????+?+??∑
∑}
1tropo iono MP 0()]L p T t t t ωδδδ??+?++ （3）
2.2 传输时间p T 计算
给定任意接收时刻，都要能够计算出传输起始时刻，从而为位置计算提供基准，这是一个模拟现实的过程[7]。本文采用迭代方法计算传输时间，更为简明实用。
首先，根据卫星星历和接收时刻计算出卫星在接收时刻的位置rec P ；
然后，由rec P 和用户的接收时刻位置user P 计算出迭代初值p T ，见式（4），c 为光速。
rec user /p T P P c =? （4）
最后，根据卫星星历和()p t T ?计算出发送起始时刻的卫星位置trans P ，再带入（4）式代替rec P 计算p T ，直到p T 收敛到要求值。
3．接收机接收信号的形成
在紧耦合导航中，为了防止失锁，由INS 对GPS 计算进行反馈，使得其能够跟踪迅速变化的频率（主要是多普勒频移变化），所以，在信号模拟中有必要对多普勒信号进行比较精细的控制：计算步长太大，使得伪距分辨率降低；计算步长太小，会使得计算变得比较繁琐，同时增大了计算量。
本文对四颗卫星信号进行仿真，采用了1ms 的计算步长，对每毫秒的数据进行计算，即每1ms 就要重新计算卫星、接收机和信道等参数，使得信号能够满足紧耦合导航的需要。之所以选用1ms 的计算步长，是因为其精度基本满足要求，同时因为伪随机码的周期为1ms ，也便于计算。
选取1ms 作为计算步长，并不是所有时候都能够满足精度要求，文献[1]中要求的计算步长为10μs 。不过，都以10μs 作为计算步长也不太合适，因为在载体运动机动性较低时，采用过小的计算步长没有必要，反而会大幅度最大计算量，此时，可以采用相对较
大的计算步长。在这里，可以使用一个选择判断，进行计算步长的自适应变化，如果机动性达到一定要求，就会采用更高精度的计算步长，否则，以一般步长进行计算。例如选取0.05m/s 的伪距变化率精度，有式（5）所表示的关系。
step step step 100m/s t 1ms
100m/s 1,000m/s t 100μs 1,000m/s 10,000m/s
t 10μs
v v v ?<=?
<<=??<<=?
径向径向径向 （5） 多颗卫星（本文中选取了四颗）经过叠加之后，再经过数字化和滤波，最后经过接收机的解调后，可以实现数据码的提取、C/A 码的捕获和跟踪，能够计算出正确的多普勒频移，从而为组合导航提供了可用的信号。尤其对于紧耦合方案，提供了可用的高动态信号。
4．仿真
4.1仿真流程
在软件信号模拟方面，实时性 关系到系统的应用。如果无法实时 完成信号的仿真，就不能使用转台 上的INS 系统，因为无法做到实时 修正GPS ，GPS 也无法实时修正 INS 。这时，可以通过对INS 系统 进行仿真，使整个系统离线运行， 验证组合导航中的紧耦合方案的可 用性。如果提高卫星、用户之间的 位置和伪距等模拟参数的计算，然 后对图1中的信号调制部分（阴影 模块）进行硬件化，或者采用效率 较高的C 代码来实现，这样就可能 满足实时性的要求。
因此，目前只完成了工作的第一 步：离线数据产生。下一步工作重点 集中在实时化方面。
4.2 仿真界面
为了实现可变参数的方便调节，以及 整个程序的可读性，建立了若干GUI 界 面，加强仿真处理的可视化、可交互性。
主界面如图2所示。主要功能是设 定若干环境参数、仿真时间和选择误差
类型等，并向用户提供仿真所处阶段，综合显示各仿真参数，绘制关键曲线。
图1 信号生成程序流程图
Fig1 Graph of the signal generation
4.3 仿真结果
该仿真程序通过接收机载体的运行模块计算出来的接收机位置和速度，加上一定的选星原则，获得选星结果之后，将伪距、开普勒频移和各项误差计算出来之后，就可以进行调制得到最终的仿真信号结果。如图3所示为仿真信号的一部分，其中：中频信号3.563MHz ，采样频率11.999MHz 。
由图3可知，仿真的信号与实际信号基本一致，其中(b)显示出GPS 信号的带宽约为2M （2×1.023M ），中频位于大约3.5M （3.563M ）处。另外，经过软件接收机的捕获、
跟踪及解算，得到的多普勒频移、卫星星历与仿真参数完全一致，说明了信号的正确性。
（a ）四颗卫星信号的时域图形（1ms 数据） （a ） Waves of the four satellites signal in time
domain(1ms)（b ） 四颗卫星信号的功率谱（1ms 数据） （b ） Power spectrum of the four satellites signal （1ms ）图3 仿真信号 Fig3 Simulation signal
图2 MATLAB 中的信号仿真交互主界面
Fig2 Mainframe of the signal simulation using MATLAB
经过GPS 软件接收机的捕获、跟踪等，依据所获得伪距、导航电文，解算出的导航误差如图4所示，获得的导航结果满足要求。
5．结论
本文以紧耦合导航中GPS 中频信号的仿真为主要目标，采用了解析表达式的仿真方式，利用真实环境参数获得GPS 信号的数字值，研究了GPS 信号仿真所需的条件和计算量的大小，为软件接收机的仿真提供了原始信号，为紧耦合技术研究奠定了基础。 不足之处是MATLAB 语言计算效率较低，尚不能够实时化仿真，所以不能够在线对该紧耦合导航系统进行仿真，如果需要进行实时化，则还需要进一步的研究工作，并将MATLAB 转化为效率更高的编程语言，利用更快的算法来实现实时应用。
参考文献
[1] Lei Dong ． IF GPS Signal Simulator Development and Verification ［D ］． UCGE Report No. 20184， 2003． [2] 赵军祥， 张其善， 常青，等． 高动态GPS 卫星信号模拟器关键技术分析及应用［J ］． 电讯技术，
2003， （4）： 49-54． [3] 赵昀， 张其善， 寇艳红． 基于Matlab 的数字中频GPS 信号生成模型［J ］． 信号处理， 2005， 21(3)：
300-303，248． [4] 张颖慧， 常青， 张其善，等． 高动态GPS 信号模拟器的软件实现［J ］． 无线电工程， 2003， 33(7)：
14-16， 60．
[5] 杨志专，胡修林，张蕴玉．基于Matlab /Simulink 的数字中频GPS 信号软件模拟器［J ］． 全球定位系统， 2006， 4： 10-13．
[6] 王惠南． GPS 导航原理与应用［M ］． 北京： 科学出版社， 2001． [7] 杨俊，武奇生． GPS 基本原理及其Matlab 仿真［M ］． 西安： 西安电子科技大学出版社， 2006． [8] Global positioning system standard positioning service signal specification. 2nd Edition ［EB/OL ］， https://www.doczj.com/doc/bd2793534.html,/pubs/gps/sigspec/gpssps1.pdf ，1995．
图4 软件接收机误差曲线
Fig4 Error curve of the receiver location
Simulation to Digital GPS IF Signal in MATLAB
Yang Yong，Chen Si，Wang Kedong
School of Astronautics, Beihang University, Beijing, PRC, (100083)
Abstract
Taking the closely coupled INS/GPS as the applied target, and according to the analysis for the GPS IF signal and the running conditions on the actual environment and the carrier, the paper gives parameters of the GPS signal delay, Doppler frequency shift, time difference, and combines the signals of many satellites using the analytic expression of the IF signal. Finally, the results in MATLAB show that signals meet the actual conditions based on the simulation calculation for the MATLAB, and proves the correctness of signal through the acquisition, track and arithmetic of the GPS software receiver.……
Keywords：GPS；high dynamic；close couple；IF signal；signal simulator
作者简介：
杨勇（1984-），男，河北人，北京航空航天大学硕士研究生，导航制导与控制专业，主要研究方向为GPS信号模拟器、地形匹配辅助导航；
陈偲（1983-），女，北京人，北京航空航天大学硕士研究生，导航制导与控制专业，主要研究方向为惯性与卫星组合导航；
王可东（1975-），男，安徽人，北京航空航天大学副教授，博士，主要从事地形匹配和组合导航技术研究。
现代数字信号处理仿真作业 1.仿真题3.17 仿真结果及图形： 图 1 基于FFT的自相关函数计算
图 3 周期图法和BT 法估计信号的功率谱 图 2 基于式3.1.2的自相关函数的计算
图 4 利用LD迭代对16阶AR模型的功率谱估计16阶AR模型的系数为： a1=-0.402637623107952-0.919787323662670i; a2=-0.013530139693503+0.024214641171318i; a3=-0.074241889634714-0.088834852915013i; a4=0.027881022353997-0.040734794506749i; a5=0.042128517350786+0.068932699075038i; a6=-0.0042799971761507 + 0.028686095385146i; a7=-0.048427890183189 - 0.019713457742372i; a8=0.0028768633718672 - 0.047990801912420i a9=0.023971346213842+ 0.046436389191530i; a10=0.026025963987732 + 0.046882756497113i; a11= -0.033929397784767 - 0.0053437929619510i; a12=0.0082735406293574 - 0.016133618316269i; a13=0.031893903622978 - 0.013709547028453i ; a14=0.0099274520678052 + 0.022233240051564i; a15=-0.0064643069578642 + 0.014130696335881i; a16=-0.061704614407581- 0.077423818476583i. 仿真程序（3_17）: clear all clc %% 产生噪声序列 N=32; %基于FFT的样本长度
基于Matlab的GMSK调制与解调 1.课程设计目的 （1）加深对GMSK基本理论知识的理解。 （2）培养独立开展科研的能力和编程能力。 （3）通过SIMULINK对BT=0.3的GMSK调制系统进行仿真。 2.课程设计要求 （1）观察基带信号和解调信号波形。 （2）观察已调信号频谱图。 （3）分析调制性能和BT参数的关系。 3.相关知识 3.1GMSK调制 调制原理图如图2.2，图中滤波器是高斯低通滤波器，它的输出直接对VCO 进行调制，以保持已调包络恒定和相位连续。 非归零数字序 高斯低通滤 波器频率调制器 （VCO） GMSK已 调信号 图3.1GMSK调制原理图 为了使输出频谱密集，前段滤波器必须具有以下待性： 1.窄带和尖锐的截止特性，以抑制FM调制器输入信号中的高频分量； 2.脉冲响应过冲量小，以防止FM调制器瞬时频偏过大； 3.保持滤波器输出脉冲响应曲线下的面积对应丁pi／2的相移。以使调制指数为1／2。前置滤波器以高斯型最能满足上述条件，这也是高斯滤波器最小移频键控(GMSK)的由来。
GMSK 信号数据 3.2GMSK 解调 GMSK 本是MSK 的一种，而MSK 又是是FSK 的一种，因此，GMSK 检波也可以采用FSK 检波器，即包络检波及同步检波。而GMSK 还可以采用时延检波，但每种检波器的误码率不同。 GMSK 非相干解调原理图如图2.3，图中是采用FM 鉴频器（斜率鉴频器或相位鉴频器）再加判别电路，实现GMSK 数据的解调输出。 图3.2GMSK 解调原理图 4.课程设计分析 4.1信号发生模块 因为GMSK 信号只需满足非归零数字信号即可，本设计中选用（Bernoulli Binary Generator）来产生一个二进制序列作为输入信号。 图4.1GMSK 信号产生器 该模块的参数设计这只主要包括以下几个。其中probability of a zero 设置为0.5表示产生的二进制序列中0出现的概率为0.5；Initial seed 为61表示随机数种子为61；sample time 为1/1000表示抽样时间即每个符号的持续时为0.001s。当仿真时间固定时，可以通过改变sample time 参数来改变码元个数。例如仿真时间为10s，若sample time 为1/1000，则码元个数为10000。 带通滤 波器限幅器判决器鉴频器GMSK 信号 输出
昆明理工大学信息工程与自动化学院学生实验报告 （2011—2012 学年第二学期） 课程名称：数字信号处理开课实验室：信自楼111 2012 年 5 月 31 日年级、专业、班生医学号姓 名 成绩 实验项目名称数字信号处理的matlab 实现指导教师 教 师 评语教师签名： 年月日 一．实验目的 熟练掌握matlab的基本操作。 了解数字信号处理的MATLAB实现。 二．实验设备 安装有matlab的PC机一台。 三．实验内容 .1.求信号x(n)=cos(6.3Пn/3)+cos(9.7Пn/30)+cos(15.3Пn/30),0≤n≤29的幅度频谱. 2. 用冲击响应不变法设计一个Butterworth低通数字滤波器,要求参数为: Wp=0.2Пαp=1dB Ws=0.3Пαs=15dB 3.用双线性变换法设计一个Chebyshev高通IIR滤波器,要求参数为: Wp=0.6Пαp=1dB Ws=0.4586Пαs=15dB 4.用窗函数法设计一个低通FIR滤波器,要求参数为: Wp=0.2Пαp=0.3dB Ws=0.25Пαs=50dB 5.用频率抽样法设计一个带通FIR滤波器,要求参数为: W1s=0.2П W1p=0.35П W2p=0.65П W2s=0.8П αs=60dB αp=1dB 6.根据 4 点矩形序列，( n ) = [1 1 1 1] 。做 DTFT 变换，再做 4 点 DFT 变换。然后分别补零做 8 点 DFT 及 16 点 DFT。 7.调用filter解差分方程，由系统对u(n)的响应判断稳定性 8编制程序求解下列系统的单位冲激响应和阶跃响应。 y[n]+ 0.75y[n -1]+ 0.125y[n -2] = x[n]- x[n -1] 四．实验源程序 1. n=[0:1:29]; x=cos(6.3*pi*n/30)+cos(9.7*pi*n/30)+cos(15.3*pi*n/30);
Matlab通信原理仿真 学号： 2142402 姓名：圣斌
实验一Matlab 基本语法与信号系统分析 一、实验目的： 1、掌握MATLAB的基本绘图方法； 2、实现绘制复指数信号的时域波形。 二、实验设备与软件环境： 1、实验设备：计算机 2、软件环境：MATLAB R2009a 三、实验内容： 1、MATLAB为用户提供了结果可视化功能，只要在命令行窗口输入相应的命令，结果就会用图形直接表示出来。 MATLAB程序如下： x = -pi::pi; y1 = sin(x); y2 = cos(x); %准备绘图数据 figure(1); %打开图形窗口 subplot(2,1,1); %确定第一幅图绘图窗口 plot(x,y1); %以x，y1绘图 title('plot(x,y1)'); %为第一幅图取名为’plot(x,y1)’ grid on; %为第一幅图绘制网格线 subplot(2,1,2) %确定第二幅图绘图窗口 plot(x,y2); %以x，y2绘图 xlabel('time'),ylabel('y') %第二幅图横坐标为’time’,纵坐标为’y’运行结果如下图： 2、上例中的图形使用的是默认的颜色和线型，MATLAB中提供了多种颜色和线型，并且可以绘制出脉冲图、误差条形图等多种形式图： MATLAB程序如下： x=-pi:.1:pi; y1=sin (x); y2=cos (x); figure (1); %subplot (2,1,1); plot (x,y1); title ('plot (x,y1)'); grid on %subplot (2,1,2); plot (x,y2);
实验6 数字滤波器的网络结构 一、实验目的： 1、加深对数字滤波器分类与结构的了解。 2、明确数字滤波器的基本结构及其相互间的转换方法。 3、掌握用MA TLAB 语言进行数字滤波器结构间相互转换的子函数及程序编写方法。 二、实验原理： 1、数字滤波器的分类 离散LSI 系统对信号的响应过程实际上就是对信号进行滤波的过程。因此，离散LSI 系统又称为数字滤波器。 数字滤波器从滤波功能上可以分为低通、高通、带通、带阻以及全通滤波器；根据单位脉冲响应的特性，又可以分为有限长单位脉冲响应滤波器（FIR ）和无限长单位脉冲响应滤波器（IIR ）。 一个离散LSI 系统可以用系统函数来表示： M -m -1-2-m m m=0 012m N -1-2-k -k 12k k k=1 b z b +b z +b z ++b z Y(z)b(z)H(z)=== =X(z)a(z) 1+a z +a z ++a z 1+a z ∑∑ 也可以用差分方程来表示： N M k m k=1 m=0 y(n)+a y(n-k)=b x(n-m)∑∑ 以上两个公式中，当a k 至少有一个不为0时，则在有限Z 平面上存在极点，表达的是以一个IIR 数字滤波器；当a k 全都为0时，系统不存在极点，表达的是一个FIR 数字滤波器。FIR 数字滤波器可以看成是IIR 数字滤波器的a k 全都为0时的一个特例。 IIR 数字滤波器的基本结构分为直接Ⅰ型、直接Ⅱ型、直接Ⅲ型、级联型和并联型。 FIR 数字滤波器的基本结构分为横截型（又称直接型或卷积型）、级联型、线性相位型及频率采样型等。本实验对线性相位型及频率采样型不做讨论，见实验10、12。 另外，滤波器的一种新型结构——格型结构也逐步投入应用，有全零点FIR 系统格型结构、全极点IIR 系统格型结构以及全零极点IIR 系统格型结构。 2、IIR 数字滤波器的基本结构与实现 （1）直接型与级联型、并联型的转换 例6-1 已知一个系统的传递函数为 -1-2-3 -1-2-3 8-4z +11z -2z H(z)=1-1.25z +0.75z -0.125z 将其从直接型（其信号流图如图6-1所示）转换为级联型和并联型。
基于MATLAB的数字信号发生器设计报告 摘要：数字信号发生器是基于软硬件实现的一种波形发生仪器。在工工程实践中需要检测和分析的各种复杂信号均可分解成各简单信号之和，而这些简单信号皆可由数字信号发生器模拟产生，因此它在工程分析和实验教学有着广泛的应用。MATLAB是一个数据分析和处理功能十分强大的工程实用软件，他的数据采集工具箱为实现数据的输入和输出提供了十分方便的函数和命令，在数字信号处理方面方便实用。本文介绍了使用MATLAB建立一个简单数字信号发生器的基本流程，并详细叙述了简单波形（正弦波、方波、三角波、锯齿波、白噪声）信号的具体实现方法。 关键字：MATLAB ，数字信号发生器 1概述 随着计算机软硬件技术的发展，越来越多现实物品的功能能够由计算机实现。信号发生器原本是模拟电子技术发展的产物，到后来的数字信号发生器也是通过硬件实现的，本文将给出通过计算机软件实现的数字信号发生器。 信号发生器是一种常用的信号源，广泛应用于电子技术实验、自控系统和科学研究等领域。传统的台式仪器如任意函数发生器等加工工艺复杂、价格高、仪器面板单调、数据存储、处理不方便。以Matlab
和LabVlEW 为代表的软件的出现，轻松地用虚拟仪器技术解决了这些问题。 Matlab 是一个数据分析和处理功能十分强大的工程实用软件，他的数据采集工具箱（data acquisition toolbox ）为实现数据的输入和输出提供了十分方便的函数和命令，利用这些函数和命令可以很容易地实现对外部物理世界的信号输出和输入。根据声卡输出信号的原理，采用Matlab 软件编程，可以方便地输出所需要的正弦波、三角波、方波等多种信号，有效地实现信号发生器的基本功能。 2 设计原理 要设计的数字信号有正弦信号、方波信号、三角波、锯齿波、白噪声、脉冲信号。其中，前五种波形都可以利用MATLAB 提供的函数实现，并根据输入的幅值、相位、频率等信息进行调整。脉冲信号由自己编写程序实现，并以定义的时间节点控制脉冲出现的时刻。 2.1 正弦信号的实现 正弦波信号的数学表达式如2.1， ()sin 2y A ft πφ=+ 2.1 其中：A 为幅值； f 为频率； φ为相位。 在MATLAB 中，相应的数字信号可以由下式2.2计算，
信源函数 randerr 产生比特误差样本 randint 产生均匀分布的随机整数矩阵 randsrc 根据给定的数字表产生随机矩阵 wgn 产生高斯白噪声 信号分析函数 biterr 计算比特误差数和比特误差率 eyediagram 绘制眼图 scatterplot 绘制分布图 symerr 计算符号误差数和符号误差率 信源编码 compand mu律/A律压缩/扩张 dpcmdeco DPCM（差分脉冲编码调制）解码dpcmenco DPCM编码 dpcmopt 优化DPCM参数 lloyds Lloyd法则优化量化器参数 quantiz 给出量化后的级和输出值 误差控制编码 bchpoly 给出二进制BCH码的性能参数和产生多项式convenc 产生卷积码 cyclgen 产生循环码的奇偶校验阵和生成矩阵cyclpoly 产生循环码的生成多项式 decode 分组码解码器 encode 分组码编码器 gen2par 将奇偶校验阵和生成矩阵互相转换gfweight 计算线性分组码的最小距离 hammgen 产生汉明码的奇偶校验阵和生成矩阵rsdecof 对Reed-Solomon编码的ASCII文件解码rsencof 用Reed-Solomon码对ASCII文件编码rspoly 给出Reed-Solomon码的生成多项式
syndtable 产生伴随解码表 vitdec 用Viterbi法则解卷积码 （误差控制编码的低级函数） bchdeco BCH解码器 bchenco BCH编码器 rsdeco Reed-Solomon解码器 rsdecode 用指数形式进行Reed-Solomon解码 rsenco Reed-Solomon编码器 rsencode 用指数形式进行Reed-Solomon编码 调制与解调 ademod 模拟通带解调器 ademodce 模拟基带解调器 amod 模拟通带调制器 amodce 模拟基带调制器 apkconst 绘制圆形的复合ASK-PSK星座图 ddemod 数字通带解调器 ddemodce 数字基带解调器 demodmap 解调后的模拟信号星座图反映射到数字信号dmod 数字通带调制器 dmodce 数字基带调制器 modmap 把数字信号映射到模拟信号星座图（以供调制）qaskdeco 从方形的QASK星座图反映射到数字信号qaskenco 把数字信号映射到方形的QASK星座图 专用滤波器 hank2sys 把一个Hankel矩阵转换成一个线性系统模型hilbiir 设计一个希尔伯特变换IIR滤波器 rcosflt 升余弦滤波器 rcosine 设计一个升余弦滤波器 （专用滤波器的低级函数） rcosfir 设计一个升余弦FIR滤波器 rcosiir 设计一个升余弦IIR滤波器 信道函数
实验5 抽样定理 一、实验目的： 1、了解用MA TLAB 语言进行时域、频域抽样及信号重建的方法。 2、进一步加深对时域、频域抽样定理的基本原理的理解。 3、观察信号抽样与恢复的图形，掌握采样频率的确定方法和插公式的编程方法。 二、实验原理： 1、时域抽样与信号的重建 （1）对连续信号进行采样 例5-1 已知一个连续时间信号sin sin(),1Hz 3 ππ=0001f(t)=(2f t)+6f t f ，取最高有限带宽频率f m =5f 0，分别显示原连续时间信号波形和F s >2f m 、F s =2f m 、F s <2f m 三情况下抽样信号的波形。 程序清单如下： %分别取Fs=fm ，Fs=2fm ，Fs=3fm 来研究问题 dt=0.1; f0=1; T0=1/f0; m=5*f0; Tm=1/fm; t=-2:dt:2; f=sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t); subplot(4,1,1); plot(t,f); axis([min(t),max(t),1.1*min(f),1.1*max(f)]); title('原连续信号和抽样信号'); for i=1:3; fs=i*fm;Ts=1/fs; n=-2:Ts:2; f=sin(2*pi*f0*n)+1/3*sin(6*pi*f0*n); subplot(4,1,i+1);stem(n,f,'filled'); axis([min(n),max(n),1.1*min(f),1.1*max(f)]); end 程序运行结果如图5-1所示：
原连续信号和抽样信号 图5-1 （2）连续信号和抽样信号的频谱 由理论分析可知，信号的频谱图可以很直观地反映出抽样信号能否恢复原模拟信号。因此，我们对上述三种情况下的时域信号求幅度谱，来进一步分析和验证时域抽样定理。 例5-2编程求解例5-1中连续信号及其三种抽样频率（F s>2f m、F s=2f m、F s<2f m）下的抽样信号的幅度谱。 程序清单如下： dt=0.1;f0=1;T0=1/f0;fm=5*f0;Tm=1/fm; t=-2:dt:2;N=length(t); f=sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t); wm=2*pi*fm;k=0:N-1;w1=k*wm/N; F1=f*exp(-j*t'*w1)*dt;subplot(4,1,1);plot(w1/(2*pi),abs(F1)); axis([0,max(4*fm),1.1*min(abs(F1)),1.1*max(abs(F1))]); for i=1:3; if i<=2 c=0;else c=1;end fs=(i+c)*fm;Ts=1/fs; n=-2:Ts:2;N=length(n); f=sin(2*pi*f0*n)+1/3*sin(6*pi*f0*n); wm=2*pi*fs;k=0:N-1; w=k*wm/N;F=f*exp(-j*n'*w)*Ts; subplot(4,1,i+1);plot(w/(2*pi),abs(F)); axis([0,max(4*fm),1.1*min(abs(F)),1.1*max(abs(F))]); end 程序运行结果如图5-2所示。 由图可见，当满足F s≥2f m条件时，抽样信号的频谱没有混叠现象；当不满足F s≥2f m 条件时，抽样信号的频谱发生了混叠，即图5-2的第二行F s＜2f m的频谱图，，在f m=5f0的围，频谱出现了镜像对称的部分。
通信原理课程设计 基于matlab的通信信道及眼图的仿真 作者: 摘要 由于多径效应和移动台运动等影响因素，使得移动信道对传输信号在时间、频率和角度上造成了色散，即时间色散、频率色散、角度色散等等，因此多径信道的特性对通信质量有着重要的影响，而多径信道的包络统计特性则是我们研究的焦点。根据不同无线环境，接收信号包络一般服从几种典型分布，如瑞利分布、莱斯分布等。因此我们对瑞利信道、莱斯信道进行了仿真并针对服从瑞利分布的多径信道进行模拟仿真。由于眼图是实验室中常用的一种评价基带传输系统的一种定性而方便的方法，“眼睛”的张开程度可以作为基带传输系统性能的一种度量，它不但反映串扰的大小，而且也可以反映信道噪声的影响。为此，我们在matlab上进行了仿真,加深对眼图的理解。 关键词：瑞利信道莱斯信道多径效应眼图 一、瑞利信道 在移动通信系统中，发射端和接收端都可能处于不停的运动状态之中，这种相对运动将产生多普勒频移。在多径信道中，发射端发出的信号通过多条路径到达接收端，这些路径具有不同的延迟和接收强度，它们之间的相互作用就形成了衰落。MATLAB中的多径瑞利衰落信道模块可以用于上述条件下的信道仿真。 多径瑞利衰落信道模块用于多径瑞利衰落信道的基带仿真，该模块的输入信号为复信号，可以为离散信号或基于帧结构的列向量信号。无线系统中接收机与发射机之间的相对运动将引起信号频率的多普勒频移，多普勒频移值由下式决定： 其中v是发射端与接收端的相对速度，θ是相对速度与二者连线的夹角，λ是信号的波长。
Fd的值可以在该模块的多普勒平移项中设置。由于多径信道反映了信号在多条路径中的传输，传输的信号经过不同的路径到达接收端，因此产生了不同的时间延迟。当信号沿着不同路径传输并相互干扰时，就会产生多径衰落现象。在模块的参数设置表中，Delay vector（延迟向量）项中，可以为每条传输路径设置不同的延迟。如果激活模块中的Normalize gain vector to 0 dB overall gain,则表示将所有路径接收信号之和定为0分贝。信号通过的路径的数量和Delay vector（延迟向量）或Gain vector(增益向量)的长度对应。Sample time(采样时间)项为采样周期。离散的Initial seed(初始化种子)参数用于设置随机数的产生。 1.1、Multipath Rayleigh Fading Channel(多径瑞利衰落信道)模块的主要参数 参数名称参数值 Doppler frequency(Hz) 40/60/80 Sample time 1e-6 Delay vector(s) [0 1e-6] Gain vector(dB) [0 -6] Initial seed 12345 使能 Normalize gain vector to 0 dB overall gain Bernoulli Random Binary Generator(伯努利二进制随机数产生器)的主要参数 参数名称参数值 Probability of a zero0.5 Initial seed54321
补充内容：模拟调制系统的MATLAB 仿真 1.抽样定理 为了用实验的手段对连续信号分析，需要先对信号进行抽样（时间上的离散化），把连续数据转变为离散数据分析。抽样（时间离散化）是模拟信号数字化的第一步。 Nyquist 抽样定律：要无失真地恢复出抽样前的信号，要求抽样频率要大于等于两倍基带信号带宽。 抽样定理建立了模拟信号和离散信号之间的关系，在Matlab 中对模拟信号的实验仿真都是通过先抽样，转变成离散信号，然后用该离散信号近似替代原来的模拟信号进行分析的。 【例1】用图形表示DSB 调制波形)4cos()2cos(t t y ππ= 及其包络线。 clf %%计算抽样时间间隔 fh=1;%%调制信号带宽(Hz) fs=100*fh;%%一般选取的抽样频率要远大于基带信号频率，即抽样时间间隔要尽可能短。 ts=1/fs; %%根据抽样时间间隔进行抽样，并计算出信号和包络 t=(0:ts:pi/2)';%抽样时间间隔要足够小，要满足抽样定理。 envelop=cos(2*pi*t);%%DSB 信号包络 y=cos(2*pi*t).*cos(4*pi*t);%已调信号 %画出已调信号包络线 plot(t,envelop,'r:','LineWidth',3); hold on plot(t,-envelop,'r:','LineWidth',3); %画出已调信号波形 plot(t,y,'b','LineWidth',3); axis([0,pi/2,-1,1])% hold off% xlabel('t'); %写出图例 【例2】用图形表示DSB 调制波形)6cos()2cos(t t y ππ= 及其包络线。 clf %%计算抽样时间间隔 fh=1;%%调制信号带宽(Hz) fs=100*fh;%抽样时间间隔要足够小，要满足抽样定理。 ts=1/fs; %%根据抽样时间间隔进行抽样
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 的数字信号处理基本操作 一、 实验目的：学会运用MA TLAB 表示的常用离散时间信号；学会运用MA TLAB 实现离 散时间信号的基本运算。 二、 实验仪器：电脑一台，MATLAB6.5或更高级版本软件一套。 三、 实验内容： （一） 离散时间信号在MATLAB 中的表示 离散时间信号是指在离散时刻才有定义的信号，简称离散信号，或者序列。离散序列通常用)(n x 来表示，自变量必须是整数。 离散时间信号的波形绘制在MATLAB 中一般用stem 函数。stem 函数的基本用法和plot 函数一样，它绘制的波形图的每个样本点上有一个小圆圈，默认是空心的。如果要实心，需使用参数“fill ”、“filled ”，或者参数“.”。由于MATLAB 中矩阵元素的个数有限，所以MA TLAB 只能表示一定时间范围内有限长度的序列；而对于无限序列，也只能在一定时间范围内表示出来。类似于连续时间信号，离散时间信号也有一些典型的离散时间信号。 1. 单位取样序列 单位取样序列)(n δ，也称为单位冲激序列，定义为 ) 0() 0(0 1)(≠=?? ?=n n n δ 要注意，单位冲激序列不是单位冲激函数的简单离散抽样，它在n =0处是取确定的值1。在MATLAB 中，冲激序列可以通过编写以下的impDT .m 文件来实现，即 function y=impDT(n) y=(n==0); %当参数为0时冲激为1,否则为0 调用该函数时n 必须为整数或整数向量。 【实例1-1】 利用MATLAB 的impDT 函数绘出单位冲激序列的波形图。 解：MATLAB 源程序为 >>n=-3:3; >>x=impDT(n); >>stem(n,x,'fill'),xlabel('n'),grid on >>title('单位冲激序列') >>axis([-3 3 -0.1 1.1]) 程序运行结果如图1-1所示。 图1-1 单位冲激序列
仿真作业 姓名：李亮 学号：S130101083
4.17程序 clc; clear; for i=1:500 sigma_v1=0.27; b(1)=-0.8458; b(2)=0.9458; a(1)=-(b(1)+b(2)); a(2)=b(1)*b(2); datlen=500; rand('state',sum(100*clock)); s=sqrt(sigma_v1)*randn(datlen,1); x=filter(1,[1,a],s); %% sigma_v2=0.1; u=x+sqrt(sigma_v2)*randn(datlen,1); d=filter(1,[1,-b(1)],s); %% w0=[1;0]; w=w0; M=length(w0); N=length(u); mu=0.005; for n=M:N ui=u(n:-1:n-M+1); y(n)=w'*ui; e(n)=d(n)-y(n); w=w+mu.*conj(e(n)).*ui; w1(n)=w(1); w2(n)=w(2); ee(:,i)=mean(e.^2,2); end end ep=mean(ee'); plot（ep）; xlabel('迭代次数');ylabel('MSE');title('学习曲线'); plot(w1); hold; plot(w2); 仿真结果：
步长0.015仿真结果 0.10.20.30.4 0.50.60.7迭代次数 M S E 学习曲线
步长0.025仿真结果
步长0.005仿真结果 4.18 程序 data_len = 512; %样本序列的长度 trials = 100; %随机试验的次数 A=zeros(data_len,2);EA=zeros(data_len,1); B=zeros(data_len,2);EB=zeros(data_len,1); for m = 1: trials a1 = -0.975; a2 = 0.95; sigma_v_2 =0.0731; v = sqrt(sigma_v_2) * randn(data_len, 1, trials);%产生v(n) u0 = [0 0]; num = 1; den = [1 a1 a2]; Zi = filtic(num, den, u0); %滤波器的初始条件 u = filter(num, den, v, Zi); %产生样本序列u(n) %（2）用LMS滤波器来估计w1和w2 mu1 = 0.05; mu2 = 0.005; w1 = zeros(2, data_len);
创新实践报告
报 告 题 目： 学 院 名 称： 姓 名：
基于 matlab 的通信系统仿真 信息工程学院 余盛泽
班 级 学 号： 指 导 老 师： 温 靖
二 O 一四年十月十五日
目录
一、引言........................................................................................................................ 3 二、仿真分析与测试 ................................................................................................... 4
2.1 随机信号的生成 ............................................................................................................... 4 2.2 信道编译码 ........................................................................................................................ 4 2.2.1 卷积码的原理 ........................................................................................................ 4 2.2.2 译码原理 ................................................................................................................ 5 2.3 调制与解调 ....................................................................................................................... 5 2.3.1 BPSK 的调制原理 .................................................................................................. 5 2.3.2 BPSK 解调原理 ...................................................................................................... 6 2.3.3 QPSK 调制与解调 ................................................................................................. 7 2.4 信道 .................................................................................................................................... 8
实验一 数字信号处理的Matlab 仿真 一、实验目的 1、掌握连续信号及其MA TLAB 实现方法； 2、掌握离散信号及其MA TLAB 实现方法 3、掌握离散信号的基本运算方法，以及MA TLAB 实现 4、了解离散傅里叶变换的MA TLAB 实现 5、了解IIR 数字滤波器设计 6、了解FIR 数字滤波器设计1 二、实验设备 计算机，Matlab 软件 三、实验内容 （一）、 连续信号及其MATLAB 实现 1、 单位冲击信号 ()0,0()1,0 t t t dt εεδδε-?=≠??=?>??? 例1.1：t=1/A=50时，单位脉冲序列的MA TLAB 实现程序如下： clear all; t1=-0.5:0.001:0; A=50; A1=1/A; n1=length(t1); u1=zeros(1,n1); t2=0:0.001:A1; t0=0; u2=A*stepfun(t2,t0); t3=A1:0.001:1; n3=length(t3); u3=zeros(1,n3); t=[t1 t2 t3]; u=[u1 u2 u3]; plot(t,u) axis([-0.5 1 0 A+2]) 2、 任意函数 ()()()f t f t d τδττ+∞ -∞=-? 例1.2：用MA TLAB 画出如下表达式的脉冲序列 ()0.4(2)0.8(1) 1.2() 1.5(1) 1.0(2)0.7(3)f n n n n n n n δδδδδδ=-+-+++++++
clear all; t=-2:1:3; N=length(t); x=zeros(1,N); x(1)=0.4; x(2)=0.8 x(3)=1.2; x(4)=1.5; x(5)=1.0; x(6)=0.7; stem(t,x); axis([-2.2 3.2 0 1.7]) 3、 单位阶跃函数 1,0()0,0t u t t ?≥?=?
% ch3example1A.m clear; f_p=2400; f_s=5000; R_p=3; R_s=25; % 设计要求指标 [n, fn]=buttord(f_p,f_s,R_p,R_s, 's'); % 计算阶数和截止频率 Wn=2*pi*fn; % 转换为角频率 [b,a]=butter(n, Wn, 's'); % 计算H(s) f=0:100:10000; % 计算频率点和频率范围 s=j*2*pi*f; % s=jw=j*2*pi*f H_s=polyval(b,s)./polyval(a,s); % 计算相应频率点处H(s)的值 figure(1); subplot(2,1,1); plot(f, 20*log10(abs(H_s))); % 幅频特性 axis([0 10000 -40 1]); xlabel('频率Hz');ylabel('幅度dB'); subplot(2,1,2); plot(f, angle(H_s)); % 相频特性 xlabel('频率Hz');ylabel('相角rad'); figure(2); freqs(b,a); % 也可用指令freqs直接画出H(s)的频率响应曲线。 % ch3example1B.m clear; f_p=2400; f_s=5000; R_p=3; R_s=25; % 设计要求指标 [n, fn]=ellipord(f_p,f_s,R_p,R_s,'s'); % 计算阶数和截止频率 Wn=2*pi*fn; % 转换为角频率 [b,a]=ellip(n,R_p,R_s,Wn,'s'); % 计算H(s) f=0:100:10000; % 计算频率点和频率范围 s=j*2*pi*f; % s=jw=j*2*pi*f H_s=polyval(b,s)./polyval(a,s); % 计算相应频率点处H(s)的值 figure(1); subplot(2,1,1); plot(f, 20*log10(abs(H_s))); % 幅频特性 axis([0 10000 -40 1]); xlabel('频率Hz');ylabel('幅度dB'); subplot(2,1,2); plot(f, angle(H_s)); % 相频特性 xlabel('频率Hz');ylabel('相角rad'); figure(2); freqs(b,a); % 也可用指令freqs直接画出H(s)的频率响应曲线。 % ch3example2A.m f_N=8000; % 采样率 f_p=2100; f_s=2500; R_p=3; R_s=25; % 设计要求指标 Ws=f_s/(f_N/2); Wp=f_p/(f_N/2); % 计算归一化频率 [n, Wn]=buttord(Wp,Ws,R_p,R_s); % 计算阶数和截止频率 [b,a]=butter(n, Wn); % 计算H(z) figure(1); freqz(b,a, 1000, 8000) % 作出H(z)的幅频相频图, freqz(b,a, 计算点数, 采样率)
数字信号处理大作业 班级：021231 学号： 姓名： 指导老师：吕雁
一写出奈奎斯特采样率和和信号稀疏采样的学习报告和体会 1、采样定理 在进行A/D信号的转换过程中，当采样频率fs.max大于信号中最高频 率fmax的2倍时(fs.max>2fmax)，采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的5～10倍；采样定 理又称奈奎斯特定理。 (1)在时域 频带为F的连续信号 f(t)可用一系列离散的采样值f(t1),f(t1±Δt)，f(t1±2Δt)，...来表示,只要这些采样点的时间间隔Δt≤1/2F，便可根据各 采样值完全恢复原始信号。 (2)在频域 当时间信号函数f(t)的最高频率分量为fmax时,f(t)的值可由一系列 采样间隔小于或等于1/2fo的采样值来确定,即采样点的重复频率fs ≥2fmax。 2、奈奎斯特采样频率 (1)概述 奈奎斯特采样定理：要使连续信号采样后能够不失真还原，采样频率必须 大于信号最高频率的两倍（即奈奎斯特频率）。 奈奎斯特频率（Nyquist frequency）是离散信号系统采样频率的一半，因哈里·奈奎斯特（Harry Nyquist）或奈奎斯特－香农采样定理得名。采样定理指出，只要离散系统的奈奎斯特频率高于被采样信号的最高频率或带宽，就可 以真实的还原被测信号。反之，会因为频谱混叠而不能真实还原被测信号。 采样定理指出，只要离散系统的奈奎斯特频率高于采样信号的最高频率或 带宽，就可以避免混叠现象。从理论上说，即使奈奎斯特频率恰好大于信号带宽，也足以通过信号的采样重建原信号。但是，重建信号的过程需要以一个低 通滤波器或者带通滤波器将在奈奎斯特频率之上的高频分量全部滤除，同时还 要保证原信号中频率在奈奎斯特频率以下的分量不发生畸变，而这是不可能实 现的。在实际应用中，为了保证抗混叠滤波器的性能，接近奈奎斯特频率的分 量在采样和信号重建的过程中可能会发生畸变。因此信号带宽通常会略小于奈 奎斯特频率，具体的情况要看所使用的滤波器的性能。需要注意的是，奈奎斯 特频率必须严格大于信号包含的最高频率。如果信号中包含的最高频率恰好为