test_gen_filter.m

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

% 本程序的目的在于：
% 已知一个滤波器频率传函的幅频特性，如何求取其时域滤波器的系数。
% 由于滤波器传函相位信息的缺失，满足该幅频特性的时域滤波器系数理论上有多个。
% 但如果需要使时域滤波器系数为实数，则可以唯一确定该滤波器系数。
% 滤波器实现的参考文献：《数字信号处理》，王世一，P243~P244
% 本程序在这里对滤波器的实现方法进行了仿真验证，验证结果完全符合理论。
clc;
clear all;
close all;
%% 获得所需的滤波器系数
% 需要的幅频分布的特性
% 归一化频率
F_Index = [-0.5 : 0.01 : 0.5];
% 归一化频率方差
Delta_F = 0.06;
% 第一种功率谱函数
% Hk_Abs = exp(-F_Index.^2/(4*(Delta_F^2)));
% 第二种功率谱函数
Fd = 0.2;
Hk_Abs = exp(-(F_Index-Fd).^2/(4*(Delta_F^2)));
% 第三种功率谱函数
% Hk_Abs = 2*a*b*sqrt(pi)*exp(-pi^2)+exp(-(F_Index-Fd).^2/(4*(Delta_F^2)));


Freq_Index = F_Index;
figure;
plot(Freq_Index, Hk_Abs);
grid on;
title('频率-幅度谱(希望的)');

% 因为前面设置的频率是-0.5到0.5，与FFT不符合，这里进行旋转，便于后续FFT计算
Hk_Abs = ifftshift(Hk_Abs);
[Hk] = Get_Hk_From_Hk_Abs(Hk_Abs);
% 获得的时域滤波器系数，确保其为实数
hk = ifft(Hk);

%% 验证产生滤波器系数的正确性
% 产生高斯白噪声
Noise_Org = (randn(1, 10000) + j* randn(1, 10000));
Noise_Later = conv(hk, Noise_Org);

% 1.利用周期图法计算原始数据功率谱，注意，对于periodogram这个函数返回的功率谱密度是以dB为单位的
[Pxx_Org_dB, w] = periodogram(Noise_Org);
[Pxx_Later_dB,w] = periodogram(Noise_Later);
Pxx_Org = 10.^(Pxx_Org_dB./10);
Pxx_Later = 10.^(Pxx_Later_dB./10);

% 2.利用pwelch方法求解功率谱密度(最终选择的就是这个方法)
nfft=256;
Fs = 1;
window=boxcar(100); %矩形窗
window1=hamming(100); %海明窗
window2=blackman(100); %blackman窗
noverlap=20; %数据无重叠
range='twosided'; % 求取双边功率谱密度
% 注意，pwelch函数返回的功率谱密度就是本身结果，并非以dB为单位，这点需要与periodogram进行区别
[Pxx_Org, w]=pwelch(Noise_Org,window2,noverlap,nfft,Fs,range);
[Pxx_Later, w]=pwelch(Noise_Later,window2,noverlap,nfft,Fs,range);

%% 显示，观察结果 
figure;
subplot(3,1,1);
Freq_Index = (-length(Pxx_Org)/2 : 1 : length(Pxx_Org)/2-1) * 1/length(Pxx_Org);
plot(Freq_Index, fftshift(Pxx_Org));
title('原始的功率谱');
grid on;
subplot(3,1,2);
Freq_Index = (-length(Pxx_Later)/2 : 1 : length(Pxx_Later)/2-1) * 1/length(Pxx_Later);
plot(Freq_Index, fftshift(Pxx_Later));
title('滤波后的功率谱');
grid on;
subplot(3,1,3);
Freq_Index = (-length(Pxx_Later)/2 : 1 : length(Pxx_Later)/2-1) * 1/length(Pxx_Later);
plot(Freq_Index, fftshift(Pxx_Later./Pxx_Org));
hold on;
Freq_Index = (-length(Hk_Abs)/2 : 1 : length(Hk_Abs)/2-1) * 1/length(Hk_Abs);
plot(Freq_Index, fftshift(Hk_Abs.^2), 'r');
grid on;
legend('仿真','理论');
title('仿真功率谱密度和理论功率谱密度');