main.m

（本文是在论坛已有一篇文章《SIRP法相干相关K分布雷达杂波的建模与仿真》的基础上修改而来

%	罗军辉，et al. Matlab 7.0在数字信号处理中的应用. 北京：机械工业出版社，2005:180-184.
%%%%K分布杂波的产生过程
clc;
clear all;
close all
%%
% K分布参数
K_Distribute_v = 2;           % 形状参数，另外一个参数在程序仿真中确定
Simulate_Data_Length=20000;   %雷达回波帧数,一帧表示一个重复周期
Nfft = 256;                         % 在程序中使用的FFT变换点数
%% 生成高斯白噪声
% xi = randn(1, Simulate_Data_Length);
% xq = randn(1, Simulate_Data_Length);
W1_k = randn(1, Simulate_Data_Length) + j * randn(1, Simulate_Data_Length);

%% 产生滤波器系数H1(k)
F_Index = [-0.5 : 1/Nfft : 0.5-1/Nfft];        % 归一化频率
Delta_F1 = 0.01;                         % 滤波器1对应的归一化频率方差
Fd1 = 0.3;                                   % 滤波器1对应的归一化多普勒中心频率
Hk1_Abs = exp(-(F_Index-Fd1).^2/(4*(Delta_F1^2)));     % 所需的滤波器1的幅频特性
Delta_F2 = 0.0001;                         % 滤波器2对应的归一化频率方差
Hk2_Abs = exp(-F_Index.^2./(4*(Delta_F2^2)));     % 所需的滤波器的幅频特性

% 求解对应的时域滤波器系数
% 因为前面设置的频率是-0.5到0.5，与FFT不符合，这里进行旋转，便于后续FFT计算
Hk_Abs = ifftshift(Hk1_Abs);
[Hk1] = Get_Hk_From_Hk_Abs(Hk_Abs);
Hk_Abs = ifftshift(Hk2_Abs);
[Hk2] = Get_Hk_From_Hk_Abs(Hk_Abs);
hk1 = ifft(Hk1);
hk2 = ifft(Hk2);
%% 生成复高斯谱杂波
% 第一种方式滤波器
%  [b,a]=butter(5,0.01);%窄带滤波器
%  Y_k=filter(b,a,W1_k);
% 第二种方式滤波器
Tmp = length(hk1);
Y_k = conv(hk1, W1_k);
Y_k = Y_k(Tmp : end);
% 根据产生杂波的平均功率求解参数，即确定了K分布的另一个参数
K_Distribute_Alpha = sqrt(var(Y_k) / (2*K_Distribute_v));

W2_k=randn(1,Simulate_Data_Length);
% 第一种方式滤波器
%  [b,a]=butter(5,0.01);%窄带滤波器
%  Z_k=filter(b,a,W2_k);
% 第二种方式滤波器
Tmp = length(hk2);
Z_k = conv(hk2, W2_k);
Z_k = Z_k(Tmp : end);

%% 通过解非线性方程，由Z(k)推求S(k)
%下面的程序解非线性方程
S_k = zeros(1, length(Z_k));

% 对Z(k)进行ZMNL变换，得到所需的功率调制变量S(k)
% 建立ZMNL变换对应的查找表
S_2_Table = (0 : 0.01 : 200);
Incomplete_Gamma_Result_Table = zeros(1, length(S_2_Table));
E_y2 = var(Y_k);
K_Distribute_Alpha_2 = K_Distribute_Alpha ^2;
for w1 = 1 : 1 : length(S_2_Table)
    tmp = E_y2 * S_2_Table(w1) / (K_Distribute_Alpha_2 * pi);
    Incomplete_Gamma_Result_Table(w1) = gammainc(K_Distribute_v, tmp);
end;
% 进行非线性方程的求解，由Z(k)推导S
Result_Tmp = 1/2 + 1/2 * erf(Z_k/sqrt(2));
S_2_Result = zeros(1, length(Z_k));
for w1 = 1 : 1 : length(Z_k)
    tmp = Result_Tmp(w1);
    tmp1 = abs(Incomplete_Gamma_Result_Table - tmp);
    [Value, Index] = min(tmp1);
    S_2_Result(w1) = S_2_Table(Index);
end;
S_k = sqrt(S_2_Result);

%% 最终合成所需的结果
X_k=Y_k.*S_k;
figure, subplot(2,1,1),plot(real(X_k));
title('K分布杂波时域波形--实部');
subplot(2,1,2),plot(imag(X_k));
title('K分布杂波时域波形--虚部');

%% 进行概率密度的检验
Num=100;
Maxdat=max(abs(X_k));
Mindat=min(abs(X_k));
NN=hist(abs(X_k),Num);
% 将每个数值区间的频率转换为概率
Pdf_Simulate = NN/(sum(NN)) / ((Maxdat-Mindat)/Num);
X_Index = (Mindat+(Maxdat-Mindat)/Num:(Maxdat-Mindat)/Num:Maxdat);
Pdf_Theory=2*((X_Index/(2*K_Distribute_Alpha)).^K_Distribute_v).*...
    besselk((K_Distribute_v-1),X_Index/K_Distribute_Alpha)./(K_Distribute_Alpha*gamma(K_Distribute_v));%theory_value
% 通过这两个式子可以计算全数据域的总概率
sum(Pdf_Simulate)*(Maxdat-Mindat)/Num
sum(Pdf_Theory)*(Maxdat-Mindat)/Num

figure,
plot(X_Index,Pdf_Simulate);
hold ;
plot(X_Index,Pdf_Theory,':red');
title('杂波的幅度分布');
xlabel('杂波的幅度')
ylabel('概率密度')
legend('估计值','理论值');
grid on;


%% 利用pwelch方法求解功率谱密度(最终选择的就是这个方法)

Fs = 1;
window=boxcar(100); %矩形窗
window1=hamming(100); %海明窗
window2=blackman(100); %blackman窗
noverlap=20; %数据无重叠
range='twosided'; % 求取双边功率谱密度
% 注意，pwelch函数返回的功率谱密度就是本身结果，并非以dB为单位
[PSD_Org, w]=pwelch(W1_k,window1,noverlap,Nfft,Fs,range);
[PSD_Later, w]=pwelch(X_k,window1,noverlap,Nfft,Fs,range);

Hk1_Simulation = PSD_Later ./ PSD_Org;
Hk1_Simulation = fftshift(Hk1_Simulation);
Hk1_Simulation = Hk1_Simulation ./ max(abs(Hk1_Simulation));
figure;
Freq_Index = (-length(Hk1_Simulation)/2 : 1 : length(Hk1_Simulation)/2-1) * 1/length(Hk1_Simulation);
plot(Freq_Index, Hk1_Simulation.^2, 'b');
% plot(Hk1_Simulation.^2, 'b');
hold on;
% figure;
Freq_Index = (-length(Hk1_Abs)/2 : 1 : length(Hk1_Abs)/2-1) * 1/length(Hk1_Abs);
plot(Freq_Index, Hk1_Abs.^2, 'r');
% plot(Hk1_Abs.^2, 'r');
legend('仿真得到','理论值');
grid on;