%完成瑞利分布隨機序列的產生，并添加角閃爍功率譜特性。 
clear all;
close all;
azi_num=25000;
fr=1000;
lamda0=0.05;
sigmav=1.0;
sigmaf=2*sigmav/lamda0;
rand('state',sum(100*clock));
d1=rand(1,azi_num);
rand('state',7*sum(100*clock)+3);%保證了兩個隨機序列的獨立性
d2=rand(1,azi_num);
xi=2*sqrt(-2*log(d1)).*cos(2*pi*d2);%產生服從高斯分布的隨機數
xq=2*sqrt(-2*log(d1)).*sin(2*pi*d2);%產生服從高斯分布的隨機數
coe_num=12;


Km = 0.2;%常用比例系數，取值范圍為0.15-0.25，典型值為0.20
Dimension = 20; %目標的最大尺寸，單位為（米）
Distance =4;
sita = Km*Dimension/Distance;%嚴格根據論文所提出的定義式寫出如右等式
sita2 = sita^2;
belda=1;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M=256; %均從圖四中的設置部分獲得相應的數據
freqx=1:1:M;
% powerfxcorr=exp(-freqx.^2/(2*sigmaf.^2));%其高斯功率譜，傅立葉變換形式不變

f_s =(2*belda*sita^2)./(belda.^2.+freqx.^2);
powerfxcorr = abs(ifft(f_s,M));
[filtercoeff,delta2] = levinson(powerfxcorr,25); 
xxi=filter(1,filtercoeff,xi);
xxq=filter(1,filtercoeff,xq);

% const=filter(1,filtercoeff,ones(1,1000));
% figure(2);
% hist(const,100);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xisigmac=std(xxi);%方差
ximuc=mean(xxi);%均值
yyi=(xxi-ximuc)/xisigmac;
xqsigmac=std(xxq);
xqmuc=mean(xxq);
yyq=(xxq-xqmuc)/xqsigmac;
sigmac=1.2;
yyi=sigmac*yyi;
yyq=sigmac*yyq;
ydata=yyi+j*yyq;
% figure(2)
% subplot(2,2,1),plot(real(ydata));title('經過濾波器后瑞利序列（實部）時域特性')
% subplot(2,2,2),plot(imag(ydata));title('經過濾波器后瑞利序列（虛部）時域特性')
% subplot(2,2,3),hist(real(ydata));title('經過濾波器后瑞利序列（實部）統計特性')
% subplot(2,2,4),hist(imag(ydata));title('經過濾波器后瑞利序列（虛部）統計特性')

num=100;
maxdat=max(abs(ydata));
mindat=min(abs(ydata));
NN=hist(abs(ydata),num);
xpdf1=num*NN/((sum(NN))*(maxdat-mindat));
xaxis1=mindat:(maxdat-mindat)/num:maxdat-(maxdat-mindat)/num;
th_val=(xaxis1./sigmac.^2).*exp(-xaxis1.^2./(2*sigmac.^2));
figure(3);plot(xaxis1,xpdf1);
hold,plot(xaxis1,th_val,'r');
title('隨機序列幅度統計結果');xlabel('幅度');ylabel('概率密度');legend('輸出隨機信號統計結果','標準瑞利分布')
signal=ydata;
signal=signal-mean(signal);%去掉直流分量


figure(1),M=256;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% psd_dat=pburg(real(signal),32,M,fr);%所得信號功率譜的計算
% psd_dat=psd_dat/(max(psd_dat));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fw = fft(ydata,M);
pw = fw.*conj(fw)/M;
N=M;
pxt(1:N/2-2) = pw(N/2+1:N-2);
pxt(N/2-1:N-4) = pw(3:N/2)/max(pw(3:N/2));
outputdisplay=pxt(N/2-1:N-4);

plot(outputdisplay,'b');
title('輸出隨機信號功率譜特性');xlabel('頻率(Hz)');ylabel('功率');
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
freqx=0:0.5*M;
freqx=freqx*fr/M;
standardf_s =(2*belda*sita^2)./(belda.^2.+freqx.^2);
plot(freqx(1:length(outputdisplay)),standardf_s(1:length(outputdisplay)),'r');hold off;
legend('輸出隨機信號功率譜特性','理想功率譜特性')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%時域

xisigmac1=std(xi);%方差
ximuc1=mean(xi);%均值
yi=(xi-ximuc1)/xisigmac1;
xqsigmac1=std(xq);
xqmuc1=mean(xq);
yq=(xq-xqmuc1)/xqsigmac1;
sigmac1=1.2;
yi=sigmac1*yi;
yq=sigmac1*yq;
ydata1=yi+j*yq;
% figure(2)
% subplot(2,2,1),plot(real(ydata));title('經過濾波器后瑞利序列（實部）時域特性')
% subplot(2,2,2),plot(imag(ydata));title('經過濾波器后瑞利序列（虛部）時域特性')
% subplot(2,2,3),hist(real(ydata));title('經過濾波器后瑞利序列（實部）統計特性')
% subplot(2,2,4),hist(imag(ydata));title('經過濾波器后瑞利序列（虛部）統計特性')

num=100;
maxdat=max(abs(ydata1));
mindat=min(abs(ydata1));
NN=hist(abs(ydata1),num);
xpdf2=num*NN/((sum(NN))*(maxdat-mindat));
xaxis1=mindat:(maxdat-mindat)/num:maxdat-(maxdat-mindat)/num;
th_val=(xaxis1./sigmac.^2).*exp(-xaxis1.^2./(2*sigmac.^2));
figure(4);
plot(xaxis1,xpdf2);
hold,plot(xaxis1,th_val,'r');
title('隨機序列幅度統計結果');xlabel('幅度');ylabel('概率密度');legend('輸出隨機信號統計結果','標準瑞利分布')