% 本程序的目的在于：
% 已知一個(gè)濾波器頻率傳函的幅頻特性，如何求取其時(shí)域?yàn)V波器的系數(shù)。
% 由于濾波器傳函相位信息的缺失，滿足該幅頻特性的時(shí)域?yàn)V波器系數(shù)理論上有多個(gè)。
% 但如果需要使時(shí)域?yàn)V波器系數(shù)為實(shí)數(shù)，則可以唯一確定該濾波器系數(shù)。
% 濾波器實(shí)現(xiàn)的參考文獻(xiàn)：《數(shù)字信號處理》，王世一，P243~P244
% 本程序在這里對濾波器的實(shí)現(xiàn)方法進(jìn)行了仿真驗(yàn)證，驗(yàn)證結(jié)果完全符合理論。
clc;
clear all;
close all;
%% 獲得所需的濾波器系數(shù)
% 需要的幅頻分布的特性
% 歸一化頻率
F_Index = [-0.5 : 0.01 : 0.5];
% 歸一化頻率方差
Delta_F = 0.06;
% 第一種功率譜函數(shù)
% Hk_Abs = exp(-F_Index.^2/(4*(Delta_F^2)));
% 第二種功率譜函數(shù)
Fd = 0.2;
Hk_Abs = exp(-(F_Index-Fd).^2/(4*(Delta_F^2)));
% 第三種功率譜函數(shù)
% 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('頻率-幅度譜(希望的)');

% 因?yàn)榍懊嬖O(shè)置的頻率是-0.5到0.5，與FFT不符合，這里進(jìn)行旋轉(zhuǎn)，便于后續(xù)FFT計(jì)算
Hk_Abs = ifftshift(Hk_Abs);
[Hk] = Get_Hk_From_Hk_Abs(Hk_Abs);
% 獲得的時(shí)域?yàn)V波器系數(shù)，確保其為實(shí)數(shù)
hk = ifft(Hk);

%% 驗(yàn)證產(chǎn)生濾波器系數(shù)的正確性
% 產(chǎn)生高斯白噪聲
Noise_Org = (randn(1, 10000) + j* randn(1, 10000));
Noise_Later = conv(hk, Noise_Org);

% 1.利用周期圖法計(jì)算原始數(shù)據(jù)功率譜，注意，對于periodogram這個(gè)函數(shù)返回的功率譜密度是以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方法求解功率譜密度(最終選擇的就是這個(gè)方法)
nfft=256;
Fs = 1;
window=boxcar(100); %矩形窗
window1=hamming(100); %海明窗
window2=blackman(100); %blackman窗
noverlap=20; %數(shù)據(jù)無重疊
range='twosided'; % 求取雙邊功率譜密度
% 注意，pwelch函數(shù)返回的功率譜密度就是本身結(jié)果，并非以dB為單位，這點(diǎn)需要與periodogram進(jìn)行區(qū)別
[Pxx_Org, w]=pwelch(Noise_Org,window2,noverlap,nfft,Fs,range);
[Pxx_Later, w]=pwelch(Noise_Later,window2,noverlap,nfft,Fs,range);

%% 顯示，觀察結(jié)果 
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('仿真功率譜密度和理論功率譜密度');