%相對對數正態分布雜波產生
clear all;close all;clc;
azi_num=2000;   %雷達回波幀數，一幀表示一個重復周期
fr=1000;    %脈沖重復頻率（Hz）
lamda0=0.05;   %波長
sigmav=1.0;
sigmaf=2*sigmav/lamda0;
rand('state',sum(100*clock));  %產生服從U（0，1）的隨機序列
d1=rand(1,azi_num);
rand('state',7*sum(100*clock)+3);
d2=rand(1,azi_num);
xi=1*(sqrt(-2*log(d1)).*cos(2*pi*d2));  %正交且獨立的高斯序列~N（0，1）；
coe_num=12;  %求濾波器系數，用傅里葉級數展開法
for n=0:coe_num
    coeff(n+1)=2*sigmaf*sqrt(pi)*exp(-4*sigmaf^2*pi^2*n^2/fr^2)/fr;
end
for n=1:2*coe_num+1
    if n<=coe_num+1
        b(n)=1/2*coeff(coe_num+2-n);
    else
        b(n)=1/2*coeff(n-coe_num);
    end
end
xxi=conv(b,xi);  %高斯譜雜波生成
xxi=xxi(coe_num*2+1:azi_num+coe_num*2);  %目的是去掉暫態響應
xsigmac=std(xxi);
xmuc=mean(xxi);
yyi=(xxi-xmuc)/xsigmac;
muc=10;   %中位值
sigmac=0.6;   %形狀參數
yyi=sigmac*yyi+log(muc);
xdata=exp(yyi);  %對數正態分布的雜波序列 
figure,plot(xdata);title('對數正態分布雜波時域波形');
%求概率密度函數的參數
num=100;
maxdat=max(abs(xdata));
mindat=min(abs(xdata));
NN=hist(abs(xdata),num);
xpdf1=num*NN/((sum(NN))*(maxdat-mindat));  %用直方圖估計的概率密度函數
xaxis1=mindat:(maxdat-mindat)/num:maxdat-(maxdat-mindat)/num;
th_val=lognpdf(xaxis1,log(muc),sigmac);
figure,plot(xaxis1,xpdf1);  %做出仿真結果的概率密度函數曲線
hold,plot(xaxis1,th_val,':r');   %做出倫理概率密度函數曲線
title('雜波幅度分布'); xlabel('幅度');  ylabel('概率密度');
signal=xdata;
signal=signal-mean(xdata);  %求功率普密度，先去掉直流分量
M=128;
psd_dat=pburg(real(signal),16,M,fr);
psd_dat=psd_dat/(max(psd_dat));   %歸一化
freqx=0:0.5*M;
freqx=freqx*fr/M;
figure,plot(freqx,psd_dat);
title('雜波頻譜');  xlabel('頻率（Hz）'); ylabel('功率普密度');
%做出理想高斯譜曲線
powerf=exp(-freqx.^2/(2*sigmaf.^2));
hold;plot(freqx,powerf,':r');