使用自功率谱、互功率谱估计滤波器幅频特性

这段时间终于对工程中的随机信号的一般处理方式有点头绪了，功率谱密度估计是十分重要的方式之一，仍需继续深入细化相关内容。

示例：使用自功率谱、互功率谱估计滤波器幅频特性，自己实现 & Matlab自带函数实现。

clc;clear;close all;
fs = 44100;
t = 0:1/fs:1-1/fs;
in = randn(size(t));
 
winlen = 1024;
overlap = winlen/2;
 
% 滤波
order = 1; % 滤波器阶数高后，则曲线不太平滑
fc = 1200;
[b,a] = butter(order, fc/(fs/2), 'low');
out = filter(b, a, in);
 
%% 自己实现cpsd、pwelch。整体思路与上周相同，但将功率谱密度估计方法由 相关法 改为 直接法
segments_i = buffer(in, winlen, overlap, "nodelay"); % 存放输入信号的分段交叠
segments_o = buffer(out, winlen, overlap, "nodelay"); % 存放输出信号的分段交叠
numSegments = size(segments_i, 2);  % 分段数
 
wd = hanning(winlen);
 
spectrum_Pxx = zeros(winlen/2+1,numSegments);  % 存放每段估计的自功率谱密度
spectrum_Pxy = zeros(winlen/2+1,numSegments);  % 存放每段估计的互功率谱密度
 
for i = 1:numSegments  % 分段估计
 
    % 估计输入信号自功率谱密度
    in_fft = fft(segments_i(:,i).*wd);
    N = length(in_fft);
 
    P1 = abs(in_fft/N);
    P2 = P1(1:floor(N/2)+1);
    P2(2:end-1) = 2*P2(2:end-1);
    Pxx_ = P2.^2/(N*fs); 
 
    % 估计输入、输出信号互功率谱密度
    out_fft = fft(segments_o(:,i).*wd);
    P3 = abs(out_fft/N);
    P4 = P3(1:floor(N/2)+1);
    P4(2:end-1) = 2*P4(2:end-1);
    Pxy_ = P4 .* conj(P2) / (N*fs);
 
    % 每段的功率谱密度存起来
    spectrum_Pxx(:,i) = Pxx_;
    spectrum_Pxy(:,i) = Pxy_;
 
end
 
f = fs*(0:N/2)/N;
 
% 计算均值
spectrumAvg_Pxx = mean(spectrum_Pxx, 2);
spectrumAvg_Pxy = mean(spectrum_Pxy, 2);
 
H2 = abs(spectrumAvg_Pxy) ./ abs(spectrumAvg_Pxx);
subplot(311);
plot(f,20*log10(H2));
title("幅频曲线(自己实现的自功率谱、互功率谱估计)");xlabel("频率");ylabel("db");
 
%% 自带pwelch、cpsd
[Pxx,f1] = pwelch(in, hanning(winlen),overlap,winlen,fs);
[Pxy,f2] = cpsd(in, out, hanning(winlen),overlap,winlen,fs);
 
subplot(312);
plot(f1, 20*log10(abs(Pxy ./ Pxx)));
title("幅频曲线(自带的pwelch、cpsd函数)");xlabel("频率");ylabel("db");
 
%% freqz理想幅频特性曲线
[H3, w] = freqz(b,a);
subplot(313);
 
f1 = w*fs/2/pi;
plot(f1,20*log10(abs(H3)));
title("幅频曲线(理想)");xlabel("频率");ylabel("db");

结果图：

遗留问题：

该示例中使用高阶滤波器后，幅频特性曲线仍会有较大的波动，有待在后续学习中解决。