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

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

示例：使用自功率谱、互功率谱估计滤波器幅频特性，自己实现 & 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_;endf = 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");

结果图：

遗留问题：

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