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

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:totalFrm

    xFrm = 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);
        
    else
        Pyx = X .* conj(X);
    end
    Pyx = Pyx(1:floor(fftLen/2) + 1);
end