一、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