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

一、welch方法原理

welch法功率谱密度估计核心总结为：信号分段、段间互相交叠、每段加窗后估算其功率谱密度，最后对所有段的估计结果进行平均得到该信号的功率谱密度。

---

二、MatLab代码实现

版本1：

clc;clear;close all;
fs = 44100;
t = 0:1/fs:1-1/fs;
in = randn(size(t));
 
winlen = 4096;
overlap = winlen/2;
 
load("myfir64.mat"); % 自己设计的64阶fir低通滤波器
filtercoe = myfir64;
out = filter(filtercoe, 1, in);

wd = hanning(winlen);
 
spectrum_Pxx = zeros(winlen/2+1,1);  % 累加每段估计的自功率谱密度
spectrum_Pxy = zeros(winlen/2+1,1);  % 累加每段估计的互功率谱密度

numSegments = fix((length(in) - overlap) / (winlen - overlap));

for i = 1:numSegments  % 分段估计
    % 自己分段帧移
    inSeg = in((i-1)*overlap+1:(i-1)*overlap+winlen);
    outSeg = out((i-1)*overlap+1:(i-1)*overlap+winlen);

    % 估计输入信号自功率谱密度
    in_fft = fft(inSeg.*wd');
    N = length(in_fft);
    in_fft = in_fft(1:floor(N/2)+1);
    in_fft(2:end-1) = 2*in_fft(2:end-1);
    Pxx_ = (abs(in_fft).^2)';
 
    % 估计输入、输出信号互功率谱密度
    out_fft = fft(outSeg.*wd');
    out_fft = out_fft(1:floor(N/2)+1);
    out_fft(2:end-1) = 2*out_fft(2:end-1);
    Pxy_ = abs(out_fft .* conj(in_fft))';
 
    % 累加每一段的功率谱密度
    spectrum_Pxx = Pxx_ + spectrum_Pxx;
    spectrum_Pxy = Pxy_ + spectrum_Pxy;
 
end

f = fs*(0:N/2)/N;
 
% 计算均值
spectrumAvg_Pxx = spectrum_Pxx ./ numSegments;
spectrumAvg_Pxy = spectrum_Pxy ./ numSegments;
 
H2 = spectrumAvg_Pxy ./ spectrumAvg_Pxx;
subplot(211);
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(212);
plot(f1, 20*log10(abs(Pxy ./ Pxx)));
title("幅频曲线(自带的pwelch、cpsd函数)");xlabel("频率");ylabel("db");

版本2：

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);

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

---

 版本3：

该版本用到如下方法：N点复序列求2个N点实序列的快速傅里叶变换

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);
figure(1);
subplot(311);
fx = w*fs/2/pi;
plot(fx, 20*log10(abs(Hx)));
title('理想的幅频特性曲线');
xlabel('频率（Hz）');
ylabel('幅值（dB）');

figure(2);
subplot(311);
phaseResponse = angle(Hx);
plot(fx, phaseResponse);
title("理想相频特性曲线");

%% 自己实现cpsd、pwelch
winlen = 4096;
overlap = winlen/2;
nfft = winlen;

inputLen = length(x);

segNum = fix((inputLen - overlap) / (winlen - overlap));
 
spectrum_Pxx = zeros(winlen/2+1,1);  % 累加每段估计的自功率谱密度
spectrum_Pxy = zeros(winlen/2+1,1);  % 累加每段估计的互功率谱密度
 
for segIdx = 1:segNum  % 估计每段

    xSeg = x((segIdx-1)*overlap+1:(segIdx-1)*overlap+winlen);  % 输入信号分段交叠
    ySeg = y((segIdx-1)*overlap+1:(segIdx-1)*overlap+winlen);  % 输出信号分段交叠
    
    Pxx_ = myCpsd(xSeg, [], hanning(winlen), nfft);    % 估计输入信号自功率谱密度
    Pxy_ = myCpsd(xSeg, ySeg, hanning(winlen), nfft);  % 估计输入、输出信号互功率谱密度

    % 累加每段的功率谱密度
    spectrum_Pxx = spectrum_Pxx + Pxx_;
    spectrum_Pxy = spectrum_Pxy + Pxy_;

end

f = fs*(0:nfft/2)/nfft;
 
% 计算均值
spectrumAvg_Pxx = spectrum_Pxx / segNum;
spectrumAvg_Pxy = spectrum_Pxy / segNum;

% 估计的幅频特性曲线
H2 = spectrumAvg_Pxy ./ spectrumAvg_Pxx;
figure(1);
subplot(313);
plot(f, 20*log10(abs(H2)));
title("自己估计的幅频曲线");xlabel("频率");ylabel("db");

% 估计的相频特性曲线
phase_response = atan2(imag(H2),real(H2));
figure(2);
subplot(312);
plot(f, phase_response);
title("自己估计的相频特性曲线");

%% matlab自带pwelch、cpsd
[Pxx,f1] = pwelch(x, hanning(winlen), overlap, winlen, fs);
[Pxy,f2] = cpsd(y, x, hanning(winlen), overlap, winlen, fs);

% 幅频特性曲线
figure(1);
subplot(312);
H5 = Pxy ./ Pxx;
plot(f1, 20*log10(abs(H5)));
title("幅频曲线(Matlab自带的pwelch、cpsd函数)");xlabel("频率");ylabel("db");

figure(2);
subplot(313);
phase_response2 = atan2(imag(H5),real(H5));
plot(f1, phase_response2);
title("相频曲线(Matlab自带的pwelch、cpsd函数)");xlabel("频率");ylabel("db");

%% 功率谱密度估计函数
function Pe = myCpsd(x, y, window, nfft)
    
    if ~isempty(y)  % 互功率谱密度

        z = x .* window' + 1i .* y .* window'; 
        % x、y序列构造为1个复数序列z

        Z = fft(z, nfft);

        X = zeros(1, nfft);
        Y = zeros(1, nfft);

        % 分别求出x、y的fft
        R = real(Z);
        I = imag(Z);

        X(1) = R(1);
        Y(1) = I(1);

        for k = 2:nfft
            X(k) = 1/2*(Z(k) + conj(Z(nfft+2-k)));   
            Y(k) = -1i/2*(Z(k) - conj(Z(nfft+2-k)));   
        end

        Pe = X .* conj(Y);
    else   % 自功率谱密度

        X = fft(x .* window', nfft);
        Pe = X .* conj(X);
    end
        Pe = Pe(1:floor(nfft/2)+1)';

end