哈工大雷达系统仿真实验报告

..雷达系统仿真实验报告：黄晓明学号：班级：1305203指导教师：俊好院系：电信学院实验一杂波和色噪声的产生—高斯谱相关对数正态随机序列的产生1、实验目的给定功率谱〔相关函数〕和概率分布，通过计算机产生该随机过程，并估计该过程的实际功率谱和概率分布以验证产生方法的有效性。2、实验原理1〕高斯白噪声的产生、均值：为位置参数、方差：、均方差：为比例参数。假设给定，那么。MATLAB中对应函数，调用根本函数产生标准正态分布。标准正态分布的产生方法有舍选抽样法、推广的舍选抽样法、变换法、极法、查表法等，其中变换法的优点是精度高，极法运算速度较变换法快1030％，查表法速度快。反变换法、反函数有理逼近法令，那么有式中，，，，，。用这一方法进展抽样，误差小于10-4。叠加法根据中心极限定理有：先产生组相互独立的上均匀分布随机数，令，那么当较大时。一般可取，那么变换对法(Box-Mullermethod)设相互独立，取，那么且相互独立。舍选法产生相互独立，令，假设满足，那么舍弃；否那么令2〕高斯色噪声的产生时域线性滤波法采用线性滤波法由高斯白噪声产生高斯色噪声的根本思想是：滤波器的暂态效应可以通过取假设干特定的输出序列初始值来解决。功率谱密度序列逆变换采用频域法反变换法来获得具有特定功率谱密度的随机序列，并基于此进一步产生高斯相关序列。设待产生的自相关函数为，对其采样得点自相关序列(注意，后半局部对应于原来的负延迟局部，类似于频谱中的负半轴。假设需产生实序列，那么满足偶函数对称性质;假设产生复序列，那么满足共轭对称)，其离散时间傅立叶变换及DFT分为：式中。设给定的连续过程/离散过程的功率谱为/，那么可通过产生一个独立随机相位序列的方法来实现该过程的总体(下两图中功率谱序列应该乘上一个系数)。该方法能获得准确的PSD期望值，具体实现方式有：频域法、时域法。3〕相关传递法相关传递法可以使一个随机序列的相关特性传递给另一个随机序列。只要使第一个序列具有所要求的振幅分布，第二个序列具有规定的相关特性，通过使第一个序列按第二个序列的大小次序排列就可使前者同时具有规定的概率密度函数和相关特性。此方法的根本思想在于：概率分布是随机序列值大小的总体描述而与其排列次序无关，而自相关特性不仅与随机序列值大小有关，更取决于序列值的相对位置，因此概率分布特性与自相关特性是两个截然不同、完全无关的概念，可以分别单独考虑实现。采用相关传递法模拟任意特定分布、功率谱的随机序列的具体步骤是：先用FFT方法得到有完全符合要求的功率谱的随机序列作为参考序列(具体实现过程见第四章高斯色噪声局部或文)，利用相关传递法对具有给定振幅分布的随机序列进展重排次序的改造，从而得到在振幅分布和谱特性两方面都符合要求的随机序列。相关传递法原理非常简单且容易实现，但其谱特性的近似程度无法保证，缺乏严格的理论分析，仅作为备用方法。4〕零记忆非线性变换法设对数正态分布密度函数：，对数正态分布y与正态分布的变换关系为：式中服从正态分布Z的相关函数式中，是、的相关系数。那么，对应的相关系数，那么待求高斯过程的相关系数为：当较小时，、，那么有，时上式说明当噪声功率较小时，对数正态分布与其对应的正态分布的相关系数近似相等。由可得到，即令，那么归一化传递函数的幅值为故产生相关对数正态分布的ZMNL法仿真框图如下：图3-3相关对数正态ZMNL法示意图具体步骤为：=1\*GB2⑴对给定杂波的功率谱密度进展次采样，建立功率谱密度序列，即(为采样点数)为了产生相关的时间序列，还需产生一个具有单位平均功率、独立的随机相位矢量序列。=2\*GB2⑵经过逆付里叶变换得到杂波的相关函数序列，利用是(3-1)给出的对应关系求出相关高斯随机序列的相关系数序列；=3\*GB2⑶产生独立高斯分布的随机序列，，且其功率谱密度。经过线性滤波器后，调制出所需要的相关正态随机序列，，其功率谱密度为；=4\*GB2⑷由参数，和非线性变换exp(.)，最终得到相关系数为服从对数正态分布的杂波序列。3、实验容给定对数正态PDF：QUOTE。取为LN(0,1)，那么QUOTE该分布均值QUOTE、方差QUOTE给定高斯型PSD：QUOTE式中为两个半功率点之间的宽度，常数使得，常数k为待定系数，与之相对应的自相关函数为功率约束关系：QUOTE，故QUOTE产生高斯白噪声；2〕产生高斯色噪声；产生高斯谱相关对数正态分布随机序列；对产生的上述随机序列进展功率谱估计；对产生的上述随机序列进展概率密度函数估计。4、实验结果及分析(1)产生高斯白噪声图1图2(2)产生高斯色噪声图3图4图5产生高斯谱相关对数正态分布随机序列；图6图7图8图9〔4〕对产生的上述随机序列进展功率谱估计；图10〔5〕对产生的上述随机序列进展概率密度函数估计。图11 分析：利用random函数产生了长度为1000，均值为0，方差为1的随机序列作为高斯白噪声，对其进展功率谱分析发现其功率谱不是理想的一条直线。考虑以下两个原因，一是因为该噪声长度有限，并不是无限大，因此必然会其频域特性产生影响。二是由于样本的随机性，一次仿真结果并不一定理想，故将其功率谱做屡次平均会得到相对较好的结果。高斯色噪声的功率谱在频域分布是不均匀的，因此将高斯白噪声通过一个线性滤波器，通过构造适当的系统函数改变其频域特性，最后经反变换得到高斯色噪声。由图5可见高斯色噪声与白噪声的区别，色噪声功率的分布不均匀对于〔3〕、〔4〕、〔5〕局部，首先设定适当参数产生功率谱密度。再由功率谱密度经傅里叶反变换得到杂波的自相关函数，利用自相关函数和高斯随机序列相关系数的对应关系得到相关系数序列。再经对应的变换得到系统函数Hw，将产生的随机序列通过该系统并做傅里叶反变换，由于方差为1，均值为零，因此直接将反变换的结果直接作为指数即得到对数正态分布序列。将0频处的分量去掉并归一化得到图10所示的功率谱密度，可见其呈高斯分布。毛刺现象较为严重，做屡次平均后可得到明显改善。利用hist命令得到该序列的概率密度函数，可见其先上升后下降并且下降速度逐渐减慢的趋势。故得到的是一个高斯相关对数正态分布序列。5、实验代码clccloseallclearallN=1000;a=2*(log(2))^0.5;fh=500;k=2*sqrt(log(2))*exp(2)/(sqrt(pi)*fh);f=1:1:N;%产生高斯白噪声n_white=random('norm',0,1,1,N);f_n_white=fft(n_white);Sn_white=abs(fftshift((f_n_white)).^2)/N;figure(1);plot(f,n_white);title('高斯白噪声');figure(2);plot(f,10*log10(Sn_white));title('高斯白噪声功率谱');xlabel('频率/Hz');ylabel('功率谱(dB/Hz)');%产生高斯色噪声Sf=k*exp(-(a*f./fh).^2);P_n_colored=n_white.*Sf;f_n_colored=sqrt(N.*Sf.*f_n_white);n_colored=ifft(f_n_colored);figure(3);plot(f,n_colored);title('高斯色噪声');axis([01000-12])figure(4);plot(f,abs(f_n_colored));title('高斯色噪声频谱');xlabel('频率/Hz');figure(5);plot(f,10*log10(P_n_colored));title('高斯色噪声功率谱');xlabel('频率/Hz');ylabel('功率谱(dB/Hz)');N=4096;%初始化各类参数a=2*(log(2))^0.5;%阿尔法fh=1200;%半功率带宽fs=N;k=2*sqrt(log(2))*exp(2)/(sqrt(pi)*fh);n=1:N/2;f=ones(1,N);f(1:N/2)=(n-1).*fs/N;f(N/2+1:N)=(N/2+1-n).*fs/N;sf=ones(1,N);%生成给定功率谱密度函数sf(1:N/2)=k*exp(-(a.*f(1:N/2)/fh).^2);sf(N/2+1:N)=k*exp(-(a.*f(N/2+1:N)/fh).^2);figure(6);plot(-N/2:N/2-1,fftshift(sf));title('模拟高斯谱相关对数正态分布序列功率谱密度函数');xlabel('频率/Hz');ylabel('幅度');s=ifft(sf);%对功率谱密度进展傅里叶反变换得到自相关函数figure(7);plot(abs(s));title('自相关函数')rou=log(s*(exp(1)-1)+1);%生成相关高斯随机序列的相关系数序列figure(8);plot(abs(rou));title('相关高斯随机序列的相关系数序列');su=fft(rou);Hw=sqrt(su/rou(1));%系统函数x=randn(1,N);%生成随机序列sx=fft(x);sy=Hw.*sx;%将生成的随机序列通过该系统y=ifft(sy);%因为均值为0，方差为1，故直接可得到z=exp(y);%高斯谱相关对数正态分布序列z=z-mean(z);%去除直流分量sz=abs(fft(z).^2)/N;%功率谱估计sz(1:10)=0;sz(4086:4096)=0;sz=fftshift(sz);sz=sz/max(sz);%进展归一化figure(9);plot(-N/2:N/2-1,sz);title('高斯谱相关对数正态分布序列功率谱密度')xlabel('频率/Hz')ylabel('幅度');figure(10);plot(abs(z));title('高斯谱相关对数正态分布序列')figure(11);hist(abs(z),0:0.1:45);title('高斯谱相关对数正态分布序列概率密度函数');实验二雷达回波的仿真与分析—阵列信号与噪声的产生及波束形成1、实验目的综合应用阵列信号模型及随机矢量产生方法，对数字波束形成技术进展原理性仿真。2、实验原理1〕独立随机矢量的产生随机矢量是由M个随机变量构成的M维随机变量，即：X=[X1X2…XM]T式中X1,X2,…,XM皆为随机变量，上标T表示矩阵转置。如果随机矢量中各随机变量相互独立，那么根据随机变量统计独立的定义，随机矢量的多维概率密度函数可表示成各随机变量PDF的乘积形式，即式中是随机变量的边沿概率密度函数由于各随机变量相互独立，因此构造此随机矢量的步骤为：1．先产生个相互独立的随机序列，其PDF为，；2．将同时刻的个随机变量直接组合起来形成列矢量即为所求。上述随机矢量模型非常适用于接收机部噪声占优的微波雷达情形，由于接收机噪声取决于各接收机部的电子噪声，与电路等元器件具体相关，因此多路具有一样构造的接收机，其部噪声是统计独立。2〕阵列远场回波模型目标点源辐射的电磁波等相位面是球面，但当目标与雷达相距足够远时，此球面由于半径相当大，其等相位面在雷达处可近似为平面，如下列图所示。3〕空间白噪声模型所谓空域白噪声是指各路噪声过程相互独立且功率相等，这对应于相互独立的随机矢量。对于空域白噪声背景，阵列的噪声协方差阵为单位阵，即，其中为噪声功率，为的单位阵。平面波平面波回波等相位面回波等相位面图1远场/Fraunhofer区情形4〕波束形成线性阵(无方向性阵元组成，连续阵可视为离散阵单元间距的特殊情形)的方向图可由加权矢量和阵列方向矢量的积表示。雷达接收天线应用领域的绝大多数情形：一般先假设目标方位，据此获得相应的匹配加权系数得到对应方位的输出，这样可以获得空间信号的方位谱，依此确定目标方位。式中，。实际阵列信号处理中依据不同准那么加权矢量有不同的表达形式。常规波束形成的加权系数通常由幅度加权及相位加权组成，相位加权为等线性变化，而幅度加权那么依据方向图旁瓣的要求，即。两头小中间大的渐变锥形幅度分布的旁射阵(不等幅但均大于零)比等幅分布的旁射阵给出较低的旁瓣，但付出了主瓣变宽的代价，这是因为在所有窗函数中矩形窗的主瓣最窄，这可利用Schelkunoff单位圆上的零点位置的变化作合理的解释。3、实验容建立等间距线性阵列模型，针对远场点目标窄带回波信号，产生阵列回波矢量及接收机噪声矢量，并进展数字波束形成仿真，定量或定性地分析不同因素变化对阵列方向图的影响。1〕分析阵元间距变化对波束形成的影响条件：等间距线性阵列，阵元数M=8，等幅加权，目标方位φr=90°，无噪声。容：分别仿真计算阵元间距d=λ/2、λ时的阵列方向图，定性分析产生的原因。2〕分析阵元数变化对波束形成的影响条件：等间距线性阵列，阵元间距d=λ/2，等幅加权，目标方位φr=90°，无噪声。容：分别仿真计算阵元数M=8、16时的阵列方向图，定性分析产生的原因〔可通过3dB波束宽度计算公式说明〕。3〕分析信噪比变化对波束形成的影响条件：等间距线性阵列，阵元数M=8，阵元间距d=λ/2，等幅加权，目标方位φr=90°，产生空间复白高斯噪声矢量。设阵列信号矢量元素s(m)=exp[jmψr]，其信号幅度为1，该阵元接收机附加的高斯噪声为n(m)=nmr+jnmi，其中实虚部均为独立同分布N(0,σ2)的高斯随机数，那么该通道合成信号x(m)=s(m)+n(m)，其噪比为：SNR=10lg[1/(2σ2)]=–3–10lg(σ2)(dB)容：分别仿真计算阵元SNR=0dB、10dB时的阵列方向图，定性分析产生的原因。4〕分析权系数变化对波束形成的影响〔选做〕条件：等间距线性阵列，阵元数M=8，阵元间距d=λ/2，等幅加权，目标方位φr=90°，无噪声。容：分别仿真计算等幅加权、汉宁加权、切比雪夫加权时的阵列方向图，定性分析产生的原因〔可分析不同窗函数引起的主瓣、旁瓣影响〕。5〕分析权系数变化对波束形成的影响〔选做〕条件：等间距线性阵列，阵元数M=8，阵元间距d=λ/2，等幅加权，无噪声。容：分别仿真计算目标方位φr=90°、60°时的阵列方向图，定性分析产生的原因〔可比拟方向图3dB主瓣宽度〕。4、实验结果及分析实验〔1〕、〔2〕局部图1图2图3实验〔3〕局部图4图5实验〔4〕局部图6图7实验〔5〕局部图8图9分析：对于实验〔1〕〔2〕局部，对照图1、图2、图3可以发现，当阵元间距一定时。阵元数越多，主旁瓣比越大，3dB波束宽度越窄，因为各阵元之间的信号是相参的而噪声噪声与噪声之间是不相关的，因此阵元数的增加相当于增加了有用信号的幅度。假设阵元间距和工作波长一样，那么就会不满足不出现栅瓣的条件，结果将出现栅瓣。对于实验〔3〕局部，当信噪比为0dB时，由图4可以发现回波信号出现了严重的失真，因为信号功率和噪声功率相当，提取有用信号困难。当信噪比为10dB时，有用信号的功率相对噪声功率较大，波形较为理想，如图5所示。对于实验〔4〕〔5〕局部，比照图6和图7可见，等幅加权时阵列方向图的旁瓣较大，进展汉宁加权后旁瓣得到了明显的抑制，但是主瓣宽度略有增加。比照发现回波入射角度越小主瓣宽度越小。5.实验代码：实验〔1〕、〔2〕局部clccloseallclearall%初始化参数。式中d中的波长会与分母中的波长%消掉，故仅用d与波长的比值表示阵元间距d=0.5;M=8;thit=90;thi=1:1:180;m=1:1:M;St=exp(j*2*pi*(m-1)*d*cos(thit*pi/180))'forangle=1:1:180S=exp(j*2*pi*(m-1)*d*cos(angle*pi/180))'f(angle)=S.'*St;endfigure(1)plot(thi,abs(f));xlabel('角度围');ylabel('幅度');title('d=0.5λ,M=8时的阵列方向图')%阵元数不变，阵元间距变为一倍波长d=1;St=exp(j*2*pi*(m-1)*d*cos(thit*pi/180))'forangle=1:180S=exp(j*2*pi*(m-1)*d*cos(angle*pi/180))'f(angle)=S.'*St;endfigure(2)plot(thi,abs(f));xlabel('角度围');ylabel('幅度');title('d=λ,M=8时的阵列方向图')%阵元数变为16个，阵元间距变为0.5倍波长d=0.5;M=16;m=1:1:M;St=exp(j*2*pi*(m-1)*d*cos(thit*pi/180))'fora=1:180S=exp(j*2*pi*(m-1)*d*cos(a*pi/180))'f(a)=S.'*St;endfigure(3)plot(thi,abs(f));xlabel('角度围');ylabel('幅度');title('d=0.5λ,M=16时的阵列方向图')实验〔3〕局部d=0.5;M=8;thit=90;snr=1;%信噪比化为直接用倍数表示sigma=1/(2)^0.5*10^(-snr*20);n=normrnd(0,sigma,M,1)+j*normrnd(0,sigma,M,1);m=1:1:M;St=exp(j*2*pi*(m-1)*d*cos(thit*pi/180))'x=St+n;fora=1:180S=hamming(8).*exp(j*2*pi*(m-1)*d*cos(a*pi/180))'f(a)=S.'*x;endfigure(1)plot(abs(f))xlabel('角度围');ylabel('幅度');title('d=0.5λ,M=8，snr=1时的阵列方向图')d=0.5;M=8;thit=90;snr=10;sigma=1/(2)^0.5*10^(-snr*20);n=normrnd(0,sigma,M,1)+j*normrnd(0,sigma,M,1);m=1:1:M;St=exp(j*2*pi*(m-1)*d*cos(thit*pi/180))'x=St+n;fora=1:180S=hamming(8).*exp(j*2*pi*(m-1)*d*cos(a*pi/180))'f(a)=S.'*x;endfigure(2)plot(abs(f))xlabel('角度围');ylabel('幅度');title('d=0.5λ,M=8，snr=10时的阵列方向图')实验〔4〕局部〔选做〕clccloseallclearall%初始化参数。式中d中的波长会与分母中的波长%消掉，故仅用d与波长的比值表示阵元间距d=0.5;M=8;thit=90;thi=1:1:180;m=1:1:M;St=exp(j*2*pi*(m-1)*d*cos(thit*pi/180))'forangle=1:1:180S=exp(j*2*pi*(m-1)*d*cos(angle*pi/180))'f(angle)=S.'*St;endfigure(1)plot(thi,abs(f));xlabel('角度围');ylabel('幅度');title('d=0.5λ,M=8时等幅加权的阵列方向图')d=0.5;M=8;thit=90;thi=1:1:180;m=1:1