功率谱密度估计 - welch方法的实现

本文主要是介绍功率谱密度估计 - welch方法的实现，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

因本人知识欠缺，后续再对下述展开讲述。

clc;clear;close all;
fs = 44100;
t = 0:1/fs:1-1/fs;
x = randn(size(t));load("myfir64.mat");
filtercoe = myfir64;
y = filter(filtercoe, 1, x);[Hx, w] = freqz(filtercoe, 1, fs);
fx = w*fs/2/pi;
subplot(211);
plot(fx, 20*log10(abs(Hx)));
title('理想幅频特性曲线');
xlabel('频率（Hz）');
ylabel('幅值（dB）');inputLen = length(x);winLen = 4096;
frmInc = winLen/2;
fftLen = winLen;
frmLen = winLen;win = hanning(winLen)';totalFrm = fix((inputLen-frmInc)/(winLen-frmInc));Pxx = zeros(fftLen/2 + 1,1);
Pyx = zeros(fftLen/2 + 1,1);% Pxx = zeros(fftLen,1);
% Pyx = zeros(fftLen,1);for frmIdx = 1:totalFrmxFrm = x((frmIdx - 1) * frmInc + 1 : (frmIdx - 1)*frmInc+frmLen);yFrm = y((frmIdx - 1) * frmInc + 1 : (frmIdx - 1)*frmInc+frmLen);PyxTmp = CPSD_feed(xFrm,yFrm,win,fftLen)';PxxTmp = CPSD_feed(xFrm,[],win,fftLen)';Pyx = Pyx + PyxTmp;Pxx = Pxx + PxxTmp;end% KMU = totalFrm*norm(win)^2;
Pyx = Pyx ./ totalFrm;
Pxx = Pxx ./ totalFrm;freqRes = 20*log10(abs(Pyx) ./ abs(Pxx));
f = fs * (0:fftLen/2)/fftLen;
subplot(212);
plot(f,freqRes);
title("估计");
% plot(freqRes);function Pyx = CPSD_feed(x,y,window,fftLen)xw = x .* window;X = fft(xw, fftLen);if ~isempty(y)yw = y .* window;Y = fft(yw, fftLen);Pyx = Y .* conj(X);elsePyx = X .* conj(X);endPyx = Pyx(1:floor(fftLen/2) + 1);
end