通信系统仿真原理与无线应用ppt课件.ppt

1,通信系统仿真原理 与无线应用,Principles of Communication SystemsSimulation with Wireless Applications,2,学习意义,设计和分析复杂的现代通信系统 ，仿 真往往必不可少的环节基于仿真的设计和分析方法成为整个通信产业中广泛使用的工具通晓仿真方法对通信领域的许多研究生的研究工作是一大支持学习仿真的学生在毕业后进入通信领域时，就已经掌握了业界所需的技能,3,第一部分 概 论,第一章 仿真的作用 现代通信系统的复杂性促进了仿真的广泛使用。这种复杂性源自现代通信系统的结构和系统运行时所处的环境。 具体表现：复杂的调制和脉冲成形技术差错控制和接收端的高级信号处理技术严格的同步技术蜂窝通信中的干扰、多径和阴影效应,4,仿真技术赖以发展的基础,数字计算机发展迅速功能强大面向通信系统的软件包的开发实用的计算能力（表现为微处理器形式）仿真理论领域的迅速发展-仿真工具、方法、技术论文和书籍,5,使用仿真的重要动机,仿真是深入理解系统特性的有价值的工具很方便地对要研究的系统进行多点测量很容易作参数研究很快观测到这些改变对系统性能的影响很容易产生时域波形、信号频谱、眼图、信号星座图、直方图和许多其他图形可以将仿真图形跟系统硬件产生的等效显示作比较 仿真的主要作用不在于获得数值而在于获得深入的理解,6,1.1 系统复杂性程度示例A. 易于解析处理的系统,7,A. 易于解析处理的系统,原因： 信道是AWGN信道，接收机是线性的 匹配滤波器是线性系统 判决统计量是高斯随机变量 假设数据源为无记忆 假设理想的符号同步意义： 可以作为更复杂的系统的基本模块； 开发出的仿真能容易得到验证,8,B. 需繁琐解析处理的系统,9,B. 需繁琐解析处理的系统,调制器与HPA后面的滤波器导致信号时间扩散，从而在时间上滤波后的信号不再局限在符号周期内，导致符号间干扰（ISI）。符号差错概率的计算将十分繁琐，经常需要用到仿真方法。 重要特性： 噪声注入点到统计量采集点之间的系统是线性的，即统计量为 Vk=Sk+Ik+Nk,10,C. 难以解析处理的系统,11,C. 难以解析处理的系统,难点： 接收机的噪声由两部分组成。 对上行链路噪声建模比较困难，即使上行链路是高斯的，判决统计量的概率密度函数还是很难确定。因此，其性能估计必须依赖于仿真。,12,1.2 仿真的多学科特点,影响研究通信系统仿真的领域线性系统理论通信原理数字信号处理(DSP)数值分析概率论估值理论计算机科学随机过程理论数论,13,1.3 仿真的作用,链路预算与系统级标校过程 链路预算主要是功率计算，包括发送功率、天线增 益、路径损耗、功率增益、放大器和滤波器的噪声 系数 关键元件的实现与测试完成硬件原型与验证仿真模型生命终结预测,14,第一部分 概 论,第二章 仿真方法论2.1 概述仿真问题涉及的步骤：将给定问题映射为仿真模型。把整个问题分解为一组小一些的问题。选择一套合适的建模、仿真和估计方法，并将其用于解决这些子问题。综合各子问题的解决结果以提供对整个问题的解决方案。,15,通信系统的基本目的是处理波形和符号 通信系统仿真试图模拟这个过程（通过产生和处理这些波形的采样值） 运行仿真的过程：用合适的输入波形驱动模型以产生输出波形分析波形以优化设计参数或获得性能指标,16,2.2 方法论的各方面,方法论中最难的步骤之一是将设计或性能估 计问题映射为仿真模型，其实现的好坏决定 了运行仿真所需的机时和仿真结果的精度。 用来解决设计或性能估计问题的整体方法或 方法论取决于具体问题的性质。 A.将问题映射到仿真模型 例：均衡器设计 系统性能评估 最终的仿真模型可通过简化初始方框图而产生。 递阶表示 划分与条件化 简化,17,B.单个模块的建模低通等效表示采样线性与非线性时不变性记忆性时域或频域仿真块处理变步长处理参数确定与其他模块的接口,18,C.随机过程建模与仿真高斯近似等效过程表示慢过程和快过程,19,2.3 性能估计,20,2.3 性能估计,BER是使用蒙特卡洛方法来确定的。BER为：差错概率的估计为：,21,第二部分 基本概念与方法,第三章 采样与量化3.1 采样低通采样定理 如果采样频率fS大于2fh，那么带限信号就可以无差错地通过其采样信号恢复。 确定性信号,22,3.1 采样,随机信号 注：为了避免混叠现象，随机信号的采样频率 仍然需要信号最高频率的2倍以上。,23,3.1 采样,带通采样定理如果带通信号的带宽为，最高频率为fh， 那么可以用大小为ffh/m的采样频率来采样并恢复信号，其中是不超过fh/的最大整数。更高的采样频率未必全部能用，除非它高于fh 。注：对于f远大于的情况，带通采样定理规 定的采样频率近似等于下界。,24,同相/正交信号采样 带通信号可表示为 注： 和 都是低通信号 和 也是低通信号,25,带通信号对应的复包络定义如下对应的频域表示为 对于xd(t) 和xq(t)的最高频率是B/2，其每个的最小采样 频率都是B。由于必须采样两个低通信号xd(t), xq(t)， 因此必须使用超过2B的采样率。,26,例题,考虑用一组采样点来表示1秒钟的某调频(FM)信号，已知载波频率为100MHz,调制信号中出现的最高频率为15KHz，试计算采样点数。（调频指数为mf =5）解：根据卡森准则可得已调信号带宽 B=2(mf +1)W=2(5+1)15=180KHz 1. 已调信号的最高频率fm=100,090KHz，根据低通采样 定理，采样1秒，需要200,180,000个采样点。 2. 现采用同相/正交形式的信号表示FM信号，其同相和正交分量的带宽为B/2，即90KHz，因而采样1秒, xd(t)和xq(t)每个至少要180,000个采样点,这样,合计1秒共有360,000个采样点。节约了200,180,000/360,000556,27,3.2 量化,设有n个量化级，每个量化级用一个长度为b比特的二进制字表示，则有 n=2b 量化信号 xqk=xk+eqk 量化信噪比 (SNR)q=S/Nq=Ex2k/Ee2qk,28,量化误差的pdf是计算机所使用的数字表示格式的函数。,定点运算 阐明量化误差的基本产生基理；运算速度比浮点计算快；定点处理器的功耗比较低。 (SNR)q=4.7712+20log10Fc+6.0206b dB浮点运算 浮点数的格式为M*(10E)，其中M和E分别是尾 数和指数 ANSI/IEEE标准规定： 尾数用53bit，指数用11bit表示,29,MATLAB 确定计算机是否采用 IEEE标准的方法,由浮点格式产生的重要参数分辨率（eps)最大数（realmax)最小数（realmin),30,% File: mparameters.m format long % 显示双精度 %a=The value of isieee is ,num2str(isieee),.; b=The value of eps is ,num2str(eps,15),.; c=The value of realmax is ,num2str(realmax,15),.; d=The value of realmin is ,num2str(realmin,15),.; %disp(a) % 显示 isieee disp(b) % 显示 eps disp(c) % 显示 realmax disp(d) % 显示 realmin format short % 恢复默认精度% End script file The value of eps is 2.22044604925031e-016. The value of realmax is 1.79769313486232e+308. The value of realmin is 2.2250738585072e-308.,31, A=1-0.4-0.3-0.2-0.1A = -2.7756e-017 结论：由浮点数算术引入的误差确实很小，在 多数的应用中可以忽略。但误差不是0， 因此计算的结果通常是不精确的。注意：在开发DSP算法时必须小心，以保证有 限字长效应对算法的影响降到最小,32,3.3 重构与内插,重构：从采样序列恢复时间连续信号方法：将采样点通过具有冲激响应h(t) 的线形滤波器，即xr(t)=xs(t)h(t)采样信号：,33,重构信号：问题：h(t)=? ，在可以接受的计算负荷范 围内获得满意的结果,34,3.3.1 理想重构,若带限信号使用超过 2fh 的采样频率进采样，则可以通过带宽为 fs/2 的理想低通滤波器来重构。,35,重构滤波器的冲激响应为,36,重构信号形式：结论：若x(t)是带限，且采样速率足够 高，则 xr(t)=x(t)，即可以得到信 号的完全重构。,37,3.3.2 上采样与下采样,直序扩频系统,38,3.3.2 上采样与下采样,问题提出：直序扩频系统中会同时出现窄带和宽带信号，如何确定采样频率？解决办法：采用两个采样率在窄带到宽带的分界处提高采样频率 通过上采样后加内插来完成在宽带到窄带的分界处降低采样频率 通过抽值来完成,39,上采样和内插,采样周期：Tu=Ts/M 采样时刻：t=nTs/M采样过程：从原有采样值x(kTs)生成新采 样值x(kTs/M)采样信号：实际信号：,40,线性内插器的脉冲响应,% File: lininterp.m function h=lininterp(M) %M=8; h1=zeros(1,(M-1); for j=1:(M-1) h1(j)=j/M; end h=0,h1,1,fliplr(h1),0; %plot(-M:M,h,go) % End of function file,41,上采样过程,1. 从xk得到xuk 2. 通过xuk和hk卷积完成内插% File: upsample.mfunction out=upsample(in_seq,M)%M=3;in_seq=1 2 3 4;L=length(in_seq);out=zeros(1,(L-1)*M+1);for j=1:L out(M*(j-1)+1)=in_seq(j);end% End function file.%out,42,上采样和内插仿真实例,% File: Ex_upsample.mM=4; % 上采样因子% 线性插值器脉冲响应h=lininterp(M); t=0:10; tu=0:40; x=sin(2*pi*t/10); xu=upsample(x,M); subplot(3,1,1);stem(t,x,bo);ylabel(x)subplot(3,1,2)stem(tu,xu,mo);ylabel(xu)xi=conv(h,xu);subplot(3,1,3)stem(xi,r.);ylabel(xi),43,下采样（抽值）,通过在采样过程中对采样周期增加M倍来完成。 采样点：xd(kTd)=x(kMTs) 采样值：xdk=xkM注意：要确保下采样信号不出现混叠现象,44,3.4 仿真采样频率,开发仿真要作的基本决定 确定采样频率影响采样频率的因素 信号带宽 系统非线性 具有反馈的系统 多径信道解决策略 选择合适的采样频率，以便在混叠误差和仿真时 间之间达成一个可以接受的折中。,45,信号传输模型信号的psd数据符号ak独立单位幅度三角波脉冲成形函数,46,仿真采样频率的确定,对于单位幅度矩形脉冲每符号m个采样(fs=m/T)混叠信噪比,47,仿真采样频率的确定,48,仿真采样频率的确定,对连续频率变量 f 进行采样, 令 f=jf1, 要求f1 1/T, 若设 f1=1/(kT), 每个旁瓣(带宽为1/T)采k个点, 则 fT=j/k 。同时, 利用 fs/2 对应km/2, 将上式中积分替换为求和，可得,49,实例：矩形脉冲波形的采样频率,% File: sna.mk=50;%每个旁瓣中采样点数nsamp=50000;%总的频率采样值snr_dB=zeros(1,20);x=1:20;for m=1:20 %每个符号采样数 signal=0;noise=0; f_fold=k*m/2; for j=1:f_fold term=(sin(pi*j/k)/(pi*j/k)2; signal=signal+term; end for j=(f_fold+1):nsamp term=(sin(pi*j/k)/(pi*j/k)2; noise=noise+term; end snr_dB(m)=10*log10(signal/noise);endplot(x,snr_dB)xlabel(每符号采样数);ylabel(混叠信噪比)% End script file,50,矩形成形脉冲和三角形成形脉冲混叠信噪比性能比较,51,作业1-11. 使用每秒10个采样的频率fs对信号 进行采样，画出X(f)和Xs(f)。假定重构滤波器是带宽为 fs/2的理想低通滤波器，具有通带增益Ts=1/fs，画出重构滤波器的输出。2. 采用 模型进行数据传送，成 形脉冲为三角形脉冲，单位幅度，周期为T。试开发 MATLAB程序画出每符号从4到20采样点时的混叠信噪比(SNR)a。和矩形成形脉冲得到的结果画在同一个坐标系里作比较，并解释其结果。,52,作业1-2QPSK信号的能量谱定义如下 其中 Tb是比特时间。试开发MATLAB程序画出每符号从4到20采样点时的混叠信噪比(SNR)a。和矩形成形脉冲得到的结果画在同一个坐标系里作比较，并解释其结果。MSK(最小移频键控)信号的能量谱定义如下 重复上面问题。,53,第四章 带通信号与系统的 低通仿真模型,4.1 带通信号的低通复包络4.1.1 复包络：时域一般的带通信号形式根据欧拉公式,54,式中 就是实信号x(t)的复包络，还可以表示成利用前面的式子，时域信号x(t)还可写成说明：通常选择f0为带通信号的中心频率， f0是 可以任意选的，可以根据方便性原则来选择。 xd(t)和xq(t) 取决于f0 ，故仿真时要选取合适的f0以便使计算量减到最小。,55,已知带通信号,56,其中Ac和fc分别表示未调载波的幅度和频率，m(t)是调制信号。假设参考时刻t0=0并假设初始相位为0。试写出该调频信号的同相和正交分量表达式。,思考题：一个模拟FM调制信号定义如下：,57,FM调制、连续时间带通模型,FM调制、连续时间低通模型,58,数字调制器的仿真模型（带通模型）,相应的离散时间信号为,59,星座图,一个信号空间被定义为由K个正交归一的基函数 所产生的K维空间，而在这个空间的信号以空间中的向量表示。对于M进制系统，在第k个信号间隔内发射的信号可表示为第m个信号的向量表示 Xm=x1m x2m , , xKm | Xm |2代表信号的能量点Xm, m=1,2,M 的K维图定义为信号的星座,60,散点图,与信号相应的散点图是xq(t)相对于xd(t) 的平面图，可按照实带通信号的低通复 包络来定义 若 xq(t)=0 或 xd(t)=0，散点图的维数为1 K=2时，散点图和星座图是紧密相关的 K2时，这种关系就消失了,61,举例：QPSK仿真结果,62,举例：16QAM仿真结果,63,4.1.2 复包络：频域,由下式可以得到两边乘以 ，经整理后得由于 是一个低通信号，故抽取上式中的低通部分,64,的傅氏变换为如果用传递函数为 的滤波器来实现低通滤波器，上式则可表示为Xd(f)与Xq(f)的表达式对Xd(f)和Xq(f)作傅氏反变换即可得到xd(t)和xq(t),65,从X(f)得到复包络频谱,66,带通信号x(t)的平均功率为低通复包络信号的平均功率为二者关系为,同相和正交分量的频谱,67,几点说明,如果实带通信号的频谱关于f0不对称，低 通复包络取复数 如果实带通信号的频谱关于f0共轭对称， 则低通复包络频谱关于f=0共轭 对称，其 包络为实数 低通复包络信号的实部和虚部都具有B/2的 带宽，为实信号带宽的一半 随机带通信号的低通复包络具有的功率是 对应的实带通信号功率的两倍 实带通信号和相应的低通复包络的信噪比 相等,68,4.2 线性带通系统,4.2.1 线性时不变系统假设系统是线性的，输入x(t)，则输出y(t)为输出的复包络表示形式为定义输入、输出的复包络满足单位冲激相应,69,根据h(t)可以确定低通复包络的频域表示即单位增益带通滤波器映射到了单位增益低通滤波器计算线性时不变系统输出的方法 用带通信号进行系统分析用复包络信号进行系统分析,70,线形带通系统的模型,71,思考题：试确定带通移相器冲激相应h(t)的复包络的同相及正交分量表达式。已知系统的输入为x(t）=Acos(2f0t+)，移相器的输出是y(t）=Acos(2f0t+)。（该模型可以用来表示解调器中同步错误）,72,4.2.2 从H(f)推导hd(t)和hq(t),由H(f)得到Hd(f)和Hq(f) ，然后对Hd(f)和 Hq(f)作 逆变换来确定hd(t)和hq(t)。 由H(f)得到 ，然后求 的逆变换来确 定 ，最后根据 确定hd(t)和hq(t)。 Hd(f)与Hq(f)的表达式,73,讨论： 如果H(f)关于f=f0共轭对称， Hq(f)和hq(t)均为 0，则 是实函数。这样滤波运算的运算负荷就可减半。 如果滤波器的带宽相对于滤波器的中心频率 较小，大部分的滤波器设计很近似地有这个 特性，则有 , 因此 可以忽略。 正交分量可以看成是H(f)关于f0的非共轭对称 性的一个度量,74,4.3 多载波信号,M个信号频分复用(FDM)的时域表示式中,75,定义y(t)的复包络为则y(t)可以表示成FDM信号的同相和正交分量分别为,76,例题 ：4路信号FDM，感兴趣的是x2(t)，令f0=f2，则合成信号y(t)的复包络为,77,注：复包络的最小采样频率取决于将哪个带通 信号平移到了f = 0处 对于f0=f2 的情况，采样频率必须满足,78,4.4 非线性于时变系统4.4.1 非线性系统,特点 非线性系统没有定义传递函数这一基本概念 非线性系统没有通过卷积将系统的输入和输出关联 非线性系统叠加不再成立方法 在对物理系统进行测量的基础上建模 通过分析来建立非线性系统的仿真模型,79,非线性系统的仿真模型（带通硬限幅器）,系统输入输入的复包络硬限幅器输出输出的复包络同项和正交分量 上式定义的器件称为带通硬限幅器，能够保证过零点位置不变的同时，消除包络上的所有变化。,80,4.4.2 时变系统,系统是线性的 时域上：通过卷积将系统的输入和输出联系起来 频域上：通过传递函数将系统的输入和输出联系起来系统是线性的、时变的 时变冲激响应 定义为在时刻t测量到的系统对此前秒在其输入端所加冲激信号的响应时变系统的传递函数,81,时变线性系统的仿真模型（两线模型）,接收信号可定义为,82,设输入信号为输入信号的复包络接收信号为,83,复路径衰减定义为接收信号可重新表示为接收信号的复包络上式定义了复低通信道模型,84,作业2-1如果信号x(t)的表达式为 （1）绘出信号及其幅度频谱的曲线。 （2）当f0=200Hz，求出其低通等效信号并绘出其幅度频谱、低 通 等效信号的同相、正交分量、包络和相位。 （3）当f0=100Hz，重做上述问题。某角度调制信号定义如下 （1）解析地确定并画出xd(t)和 xq(t) 。 （2）使用MATLAB和FFT变换，确定并画出X(f)，要求画出幅度 和相位。 （3）使用MATLAB，确定并画出等效低通信号的频谱。 （4）使用MATLAB，确定并画出xd(t)和 xq(t) 。,85,作业2-2（选作）3. 一个带通硬限幅器的输入为AM信号 硬限幅器的输出为 （1）使用MATLAB，确定x(t)的复包络 ，并画出xd(t)和xq(t) 。 （2）编写MATLAB仿真程序来实现带通硬限幅器的仿真模型。 （3）使用 作为仿真模型的输入，产生并画出yd(t)， yq(t) 和 。 （4）使用仿真产生的yd(t)， yq(t) ，产生并画出y(t)和，将所得的 结果和y(t)解析表达式作比较 。 （5）求解xd(t)，xq(t)， yd(t)和 yq(t)的解析表达式 ，画出对应的曲 线，并和利用仿真方法得到的结果作比较。,86,第五章 随机信号的产生与处理,研究目的：随机数发生器在通信系统仿真中的运用 研究内容：产生和测试描述随机过程的随机变量序列 伪随机序列：（pseudo-random sequence）在观测（仿真）区间上“呈现随机性”的序列 噪声波形模型精度：取决于具体的应用 研究工具：MATLAB,87,5.1 平稳与遍历过程,在仿真通信系统时所作的前提假设条件： 用于表示信号、噪声和干扰而产生的样本函 数为各态历经的 仿真所处理的波形是由其内在的统计模型定 义的总体（Ensemble）中的一个典型成员 仿真计算得到的时间平均等于总体平均 上述假设条件对应于随机过程具有遍历性。 遍历过程总是平稳的，因此可以假设仿真中 产生的样本函数属于平稳随机过程。,88,举例说明 三组样本函数,89,观察结果,1. Ex(t)=0；当i时，时间平均收敛于02. 当t=0.875和1.875时， Ey(t) 1 当t=0.375和1.375时， Ey(t)-1 当t=0.125、0.625、1.125和1.625时，Ey(t)03. 若以t=0.5k对z(t)进行采样，当k为偶数时， Ez(t)1.5，当k为奇数时， Ez(t)-1.5 。4. y(t) 和z(t) 为周期平稳过程，即其矩是周期的,90,数字调制器中一个常用模块：random_binary,% File: function random_binary function x,bits=random_binary(nbits,nsamples) x=zeros(1,nbits*nsamples) bits=round(rand(1,nbits); for m=1:nbits for n=1:nsamples index=(m-1)*nsamples+n; x(1,index)=(-1)bits(m); end end % End of function 例：QPSK调制 x= random_binary(nbits,nsamples)+j* random_binary(nbits,nsamples),91,应用实例：产生一个10比特的QPSK信号，采样频率为每比特8个采样点,92,5.2 均匀随机数发生器5.2.1 线性同余,线性同余发生器 (linear congruence generator, LCG) xi+1=axi+cmod(m) x0 称为种子常用算法混合同余法具有素模数的乘性算法具有非素模数的乘性算法,93,混合同余法 xi+1=axi+cmod(m) , c0时，最大周期为m满足条件： 增量c与m互质，即c与m没有素公因子 a-1是p的倍数，p为m的任意一个素因子 如果m是4的倍数，则a-1是4的倍数例题：设计一个周期m=5000的混合同余发生器1.设定c 5000=(23)(54) 令c等于2和5之外的素数的乘积 可选 c=(32)(72)=1323 2.设定a a必须满足 a-1=k1p1(p1=2)和a-1=k2p2 (p2=5) 因为4是m的因数，则有a-1=4k3，显然可令 a-1=245k=40k，可取k=6，则a=2413.结果 xi+1=241xi+1323mod(5000) 全周期发生器,94,5.2.2 随机数发生器的测试-利用散点图 xi+1=(65xi+1mod(2048) xi+1=1229xi+1mod(2048),a=65 a=1229,95,5.2.2 随机数发生器的测试-Durbin-Watson xi+1=(65xi+1mod(2048) xi+1=1229xi+1mod(2048),Durbin-Watson独立性测试可以通过计算Durbin参数来完成 讨论：D=2(1-)，即有0D4。 =0时，D=2 表示不相关 D2 表示负相关 结果：当a= 65时， D=1.9925, =0.0037273 当a=1229时， D=1.6037, =0.19814,96,5.2.3 最低标准表明给定LCG能通过各种随机性的统计测试，以便彻底解决测试它的质量。,全周期能通过所有适用的随机性统计测试易于从一台计算机移植到另一台计算机Lewis,Goodman和Miller最低标准 xi+1=(16807xi)mod(2147483647)最低标准均匀随机数发生器 Wichmann-Hill算法,97,5.2.4 MATLAB实现,函数库内均匀随机数发生器 rand 种子数与种子向量（35个元素） 注意事项可以使用默认的种子数，也可以自己定义种子数可以产生相同的随机数序列种子数可以默认，也可自己给定系统时钟可以用来随机化初始的种子数种子数存在一个缓冲区，而不是工作区，因此，clear all 对种子数不起作用,98,5.3 将均匀分布随机变量映射成任意pdf,已知目标随机变量的概率分布函数(CDF)具有闭合形式 采用逆变换法已知目标随机变量的pdf具有闭合形式，但CDF不具有闭合形式 采用舍弃法已知pdf和CDF都不具有闭合形式。 采用直方图法,99,5.3.1 逆变换法,可以将一个不相关的均匀分布的随机序列U变换为一个具有概率分布函数Fx(x)的不相关的(独立的样本)序列X。使用此法要求已知分布函数Fx(x)为闭合形式。例1 将一均匀分布的随机变量变换成具有单边指数分布 的随机变量 , 其中u(x) 是单位阶跃函数。 随机变量U 与1-U 等价,100,逆变换法概率分布函数,101,举例：均匀分布转换为指数分布,102,举例：均匀分布转换为指数分布,103,例2 将一均匀分布的随机变量变换成具有瑞利分布的随机变量 , 其中u(x)是单位阶跃函数。随机变量U 与1-U 等价,104,举例：均匀分布转换为瑞利分布,105,5.3.2 直方图法,产生实验数据的直方图 根据直方图大体了解pdf和CDF 运用逆变换法求得随机变量原理：样本值落在直方图第i个直方上的概率为在x点处估算的CDF值可表示为式中令FX(X)=U，其中U是均匀随机变量，代入上式，解得X,106,直方图法实现步骤,从一个能产生(0,1)区间均匀分布的随机 数发生器中取出样本，产生随机变量U 确定满足如下条件的i值 Fi-1UFi 根据X公式产生随机变量序列X,107,5.3.3 舍弃法思想,用函数MgX(x)界定目标pdf，其中可设定gX(x)为 均匀分布在(0,a)区间，M为常数且满足 MgX(x) fX(x),108,图示,109,算法步骤,1. 产生在(0, 1)上均匀分布的U1和U2。2. 产生在(0, a)上均匀分布的V1，其中a是X 的最大值。3. 产生在(0, b)上均匀分布的V2，其中b是不 小于fx(x)的最大值4. 如果V2 fx(V1) ，令X=V1；如果不等式不满 足，则丢弃V1和V2 ，从第一步开始重复以 上过程 通过以上算法可以产生具有目标pdf fx(x)的随机变量,110,例3 利用舍弃法产生具有如下pdf的随机变量X 令a=R，b=M/R N=10000 接受样本数=7809 舍弃样本数=2191 效率=78.1%,111,112,5.4 产生不相关的高斯随机变量,高斯随机变量的CDF其中Q(x)定义为注意：由于高斯Q(x)不能写成闭合形式，所以无法使用逆变换方法，尽管可以用舍弃法，但效率不高。,113,5.4 .1 均匀变量求和法,理论依据：中心极限定理(central limit theorem)描述内容：广义条件下，当N时，N个独立随机变量之和的pdf收敛为高斯随机变量的pdf。构成方法：,114,讨论： 由Ui-1/2的变化范围是(-1/2, 1/2)，可知Y的值 域为(-BN/2, BN/2)，从而导致pdf的拖尾被截 短到BN/2。 当确定数字通信系统误码率时， pdf的拖尾 是十分重要的，因为pdf的拖尾表示导致传输 差错的大噪声。 由 推出 可知高斯随机变量只在上述范围取非零值。 当N=100时，随机变量在17.32y处被截短,115,数字通信系统性能和噪声pdf截短拖尾间的关系,接收器输入端低SNR接收器输入端高SNR,116,优点：在Y的均值附近能很好地近似高斯随机变量 即使组成Y 的随机变量Ui 的pdf不是均匀分布的，Y 照样可以近似为高斯的算法的选择： 不适合于传统串行处理应用的算法可能会相当适合用并行处理机来处理,117,5.4 .2 瑞利随机变量映射到高斯随机变量,联合pdf变换成极坐标,118,边界概率密度R是瑞利RV，是在(0,2)的均匀RV构造高斯随机变量：均值为0，方差为2X=Rcos Y=Rsin ,119,5.4 .2 瑞利随机变量映射到高斯随机变量算法步骤：Box-Muller算法,1. 利用均匀随机变量U生成瑞利随机变量R2. 产生一个在(0，2)上均匀分布的随机变量3. 利用瑞利随机变量的正交投影产生一对高斯随机变量 X和Y等效为：一对独立的高斯随机变量可由一对均匀分布的 随机变量U1和U2通过如下算式产生,120,5.4 .3 极坐标法算法步骤：,1.产生在区间(0，1)上分布的两个独立的随机变 量U1和U2。2.令V1=2U1-1， V2=2U2-1，使得V1和V2在区间 (-1，1)上为均匀分布。3.令 。如果S1，进入第四步；若果 S1，舍弃S值并跳回第一步。4.令 。5.令X=A(S) V1，Y=A(S) V2 。注：极坐标算法得到结果的相关性要比Box-muller算法好一些；但要产生一定数目的高斯随机变量，需要调用极坐标算法的次数未知，其缺点会使仿真程序变复杂。,121,5.5 产生相关的高斯随机变量5.5.1 确定给定的相关系数,构造方法：1. 产生一对不相关的高斯随机变量X和Y2. 产生另一个随机变量Z，定义如下3. X和Z都是零均值且具有相同方差的随 机变量。可以证明X和Z的相关系数是 。,122,证明：Z的均值Z的方差XZ的协方差相关系数,123,构造复的相关随机过程算法步骤,产生两组相互独立的零均值，方差相同的Gaussian RV，即x1(t), y1(t)，x2(t), y2(t)构造两组相关系数为的RV，x1(t), z1(t)， x2(t), z2(t)构造相关系数为的复随机序列X0和 Z0 X0= x1(t)+j x2(t) Z0 =z1(t)+jz2(t)思考题：试证明随机变量X0和 Z0的相关系数为,124,5.5.2 确定任意的功率谱密度或自相关函数,思路：对一组不相关的样本进行适当的滤波，从而使之具有目标功率谱密度。不相关样本的功率谱密度在仿真带宽| f |fs/2上 是常数，方差为 随机数（噪声）的方差和采样频率的关系 对于独立样本，其功率谱密度为常数 K=N0/2Watts/Hz 为得到目标功率谱密度，所需的H(f )为,125,随机过程通过线性系统独立样本的PSD相关随机序列的产生,126,举例 1：设计具有给定传递函数的滤波器,方法：确定传递函数在最小均方误差意义下的 最佳拟合。可通过MATLAB函数yulewalk 求解Yule-Walker方程来完成。实现：利用yulewalk求出滤波器的系数bk和ak， 使得传递函数 是给定传递函数的最小均方误差拟合。,127,假定给定的功率谱密度具有S(f)=K/f的形式（闪烁噪声），用来建模振荡器的相位噪声。要产生这种噪声，必须使用一个滤波器，使其传递函数为 。因为当f0时,H(f),故传递函数定义为利用MATALAB实现H(f)，其中频率用Nyquist频率做归一化处理，即有f0=fN/20。,128,f=0:200; % frequency pointsfn=200; % Nyquist rateF=f./fn; % frequency vectorM=abs(1./sqrt(F); % normalized frequency responseM=zeros(1,11),M(11:200); size(M) % bound from zero frequencyb1,a1=yulewalk(3,F,M); % generate order=3 filterb2,a2=yulewalk(20,F,M); % generate order=20 filterh1,w1=freqz(b1,a1); % generate 3-rd order H(f)h2,w2=freqz(b2,a2); % generate 20-rd order H(f)subplot(2,1,1)plot(F,M,:,w1/pi,abs(h1)xlabel(归一化频率)ylabel(幅度响应)subplot(2,1,2)plot(F,M,:,w2/pi,abs(h2)xlabel(归一化频率)ylabel(幅度响应),129,图中虚线表示期望特性，实线表示近似特性上图为3阶近似，下图为20阶近似,130,举例 2:用于无线移动系统仿真的随机过程,所求滤波器的传递函数为,131,% Generate and test the impulse response of the filter.clear all;clc;fd=100; % maximum dopplerimpw=jakes_filter(fd); % call to Jakes filterfs=16*fd;ts=1/fs;time=1*ts:ts:128*ts;subplot(2,2,1)stem(time,impw,.);grid;xlabel(Time);ylabel(Impulse Response);% -% Square the FFT and check the power transfer functionh f=linear_fft(impw,128,ts); % generate H(f) for filtersubplot(2,2,2)plot(f,abs(h.*h);grid;xlabel(Frequency);ylabel(PSD),132,% Put Gaussian noise through and check the output psdx=randn(1,1024); % generate Gaussian inputy=filter(impw,1,x); % filter Gaussian inputoutput_psd ff=log_psd(y,1024,ts); %log of PSD figuresubplot(2,2,3)；plot(ff,output_psd);grid;axis(-500 500 -50 0)xlabel(Frequency);ylabel(PSD);% FIlter complex noise and look at the envelope fdaingz=randn(1,1024)+j*randn(1,1024); zz=filter(impw,1,z); % filter the complex noisetime=(0.0:ts:1024*ts); % new time axis% Normalize output and plot envelopezz=zz/max(abs(zz); % normalize to onesubplot(2,2,4)plot(time(161:480),10*log10(abs(zz(161:480);grid;axis(0.1 0.3 -20 0)xlabel(Time);ylabel(Log Amplitude),133,function impw=jakes_filter(fd)% FIR implementation of Ja