通过zmnl方法，完成了对瑞丽杂波，威布尔杂波，对数正态分布杂波以及K杂波的建模仿真

clear all;close all;clc; N=8000; alf=1.5; v=1; fr=1000; %f3db=40; %f=-f3db:2*f3db/N:f3db; %S=exp(-(1.665*f/f3db).^2);%构造高斯功率谱 %s=ifft(S(1:N));%求得傅里叶逆变换也就是相关系数 num=8000;%采样点数数据点数 n=0:1:num-1; %-----由功率谱和相关系数确定相关高斯序列； fs=1000;%采样频率 f=fs/num*(-num/2:num/2-1);%对称的功率谱 sigmaf=50; amp=1/(sqrt(2*pi)*sigmaf); Pw=amp*exp(-(f).^2/(2*sigmaf^2));%高斯功率谱 sn=ifft(Pw);%由最终的高斯功率分布得到相关系数序列 s1=0:0.1:1;%数据源的纵坐标 从名称命名上看应该是与上边的相关系数对应 但是我细致的看了一下后面的程序这个貌似应该是输入杂波的相关系数比较合理 首先是拟合出来的函数应该是以输出的相关系数为已知变量（横坐标）输入的为未知变量的这样一个函数p1=polyfit(y,s1,4)以y为横坐标，s1为纵坐标其次是确定数值比较方便直接带入数值即可 %y=[0 0.0034 0.0134 0.0302 0.0538 0.0844... % 0.1224 0.1682 0.2278 0.3536 1.0000];%数据源的横坐标 y=[0 -0.0005 -0.0020 -0.0045 -0.0080 -0.0126 -0.0181 -0.0248 -0.0331 -0.0477 -0.1050]; %figure(2) q=0:1:10; plot(y,s1,'*') p1=polyfit(y,s1,4);%曲线4阶拟合 得到的是拟合曲线多项式的系数 得到的位数是5正好是最高的次数4+1 x=polyval(p1,y);%检验拟合效果 plot(y,x) r=polyval(p1,sn);%多项式在s处的值 s是最终杂波的相关系数 也就是此时他所对应的其中一个输入杂波的相关系数 pp=r.^10;%求另外一个相关系数 Xw1=fft(r,num); %相关高斯随机序列功率谱 abs_Hw1=sqrt(abs(Xw1)/r(1));%这个有一个流程图可以表示 Xw2=fft(pp,num); abs_Hw2=sqrt(abs(Xw2)/pp(1)); w1=rand(1,num);%生成随机序列但是均匀分布的 w2=rand(1,num); ww1=2*pi*w1; ww2=2*pi*w2; ww1=ww1/max(abs(ww1));%应该是归一化用 ww2=ww2/max(abs(ww2)); nois1=exp(j*(ww1)); %加入的随机相位 生成第一个滤波器 Hw1=abs_Hw1.*nois1; %成形滤波器频谱 nois2=exp(j*(ww2));%生成第二个滤波器 Hw2=abs_Hw2.*nois2; vi1=randn(1,8000);%得到高斯白噪声 vi2=randn(1,8000); vi3=randn(1,8000); vi4=randn(1,8000); %vi5=randn(1,8000); %vi6=randn(1,8000); %vi7=randn(1,8000); %vi8=randn(1,8000); %vi9=randn(1,8000); %vi10=randn(1,8000); vq1=randn(1,8000);%得到高斯白噪声 vq2=randn(1,8000); vq3=randn(1,8000); vq4=randn(1,8000); %vq5=randn(1,8000); %vq6=randn(1,8000); %vq7=randn(1,8000); %vq8=randn(1,8000); %vq9=randn(1,8000); %vq10=randn(1,8000); v1=vi1+j*vq1; v2=vi2+j*vq2; v3=vi3+j*vq3; v4=vi4+j*vq4; subplot(2,2,1) plot(abs(v1)) title('独立不相关高斯随机序列') xlabel('a') %v5=vi5+j*vq5; %v6=vi6+j*vq6; %v7=vi7+j*vq7; %v8=vi8+j*vq8; %v9=vi9+j*vq9; %v10=vi10+j*vq10; %figure(3) %plot(p,v1) %title('独立不相关高斯随机序列') %y1=xcorr(v1);%求解自相关 %Z1=fft(y1(1:num+1));%自相关序列进行傅里叶变换应该就是求4解z的功率谱 %figure(5) %plot(Z);%做出来的图像不是很理解 %Z1=10*log10(abs(Z1(1:num/2+1))/sum(abs(Z1)));%这里是在做什么 求解什么参数？ %figure(5) %plot(Z1); %v1=(v1-mean(v1))/std(v1)*alf^2; %v2=(v2-mean(v2))/std(v2)*alf^2; %v3=(v3-mean(v3))/std(v3)*alf^2; %v4=(v4-mean(v4))/std(v4)*alf^2; vi11=randn(1,8000); vq11=randn(1,8000); vi12=randn(1,8000); vq12=randn(1,8000); v11=vi11+j*vq11; v12=vi12+j*vq12; %v11=(v11-mean(v11))/std(v11); %v12=(v12-mean(v12))/std(v12); V1=fft(v1); V2=fft(v2); V3=fft(v3); V4=fft(v4); %V5=fft(v5); %V6=fft(v6); %V7=fft(v7); %V8=fft(v8); %V9=fft(v9); %V10=fft(v10); V11=fft(v11); V12=fft(v12); W1=Hw1.*V1;W2=Hw1.*V2;W3=Hw1.*V3;W4=Hw1.*V4;%W5=Hw1.*V5; %W6=Hw1.*V6;W7=Hw1.*V7;W8=Hw1.*V8;W9=Hw1.*V9;W10=Hw1.*V10; W11=Hw2.*V11;W12=Hw2.*V12; w1=ifft(W1);w2=ifft(W2);w3=ifft(W3);w4=ifft(W4);%w5=ifft(W5); %w6=ifft(W6);w7=ifft(W7);w8=ifft(W8);w9=ifft(W9);w10=ifft(W10); w11=ifft(W11);w12=ifft(W12); std(w12) w1=(w1-mean(w1))/std(w1)*alf^2; w2=(w2-mean(w2))/std(w2)*alf^2; w3=(w3-mean(w3))/std(w3)*alf^2; w4=(w4-mean(w4))/std(w4)*alf^2; std(w1) mean(w1) w11=(w11-mean(w11))/std(w11); w12=(w12-mean(w12))/std(w12); std(w12) mean(w12) w1=w1.^2;w2=w2.^2;w3=w3.^2;w4=w4.^2;%w5=w5.^2; %w6=w6.^2;w7=w7.^2;w8=w8.^2;w9=w9.^2;w10=w10.^2; w11=w11.^2;w12=w12.^2; Y=w1+w2+w3+w4; %Y=w1+w2+w3+w4+w5+w6+w7+w8+w9+w10; X=w11+w12; ydata=X+Y; ydata=sqrt(ydata); subplot(2,2,2) plot(abs(ydata)) title('K分布杂波') xlabel('b') %--------------------- num_pdf=100; maxdat=max(abs(ydata)); mindat=min(abs(ydata)); NN=hist(abs(ydata),num_pdf);%绘制直方图 xpdf1=num_pdf*NN/((sum(NN))*(maxdat-mindat));%幅度分布 xaxis1=mindat:(maxdat-mindat)/num_pdf:maxdat-(maxdat-mindat)/num_pdf;%横坐标 th_va1=2./(alf*gamma(v)).*((xaxis1/(2*alf)).^(v)).*besselk(v-1,xaxis1/alf);%理论幅度分布 subplot(2,2,3),plot(xaxis1,xpdf1); %作出仿真结果的概率密度函数曲线 hold ,plot(xaxis1,th_va1,':r'); %作出理论概率密度函数曲线 title('杂波的幅度分布'); xlabel('杂波的幅度') ylabel('概率密度') %%%%用burg法来估计功率谱密度%%%%% signal=ydata; signal=signal-mean(signal); %求功率谱密度，先去掉直流分量 到中心 M=256; psd_dat=pburg(signal,32,M,fr); psd_dat=psd_dat/(max(psd_dat)); %归一化处理 psd_dat=[psd_dat(126:256)',psd_dat(1:125)']'; freqx=-0.5*M:1:0.5*M-1; freqx=freqx*fr/M; subplot(2,2,4); plot(freqx,psd_dat); title('杂波频谱'); xlabel('频率/Hz') ylabel('功率谱密度') %%%作出理想的功率谱曲线%%% powerf=exp(-freqx.^2/(2*sigmaf.^2)); hold plot(freqx,powerf,':r');