pwelch
韦尔奇的功率谱密度估计
语法
说明
pxx = pwelch(x)x 的功率谱密度 (PSD) 估计值 pxx。如果 x 是向量，则它被视为单通道。如果 x 是矩阵，则会针对每列独立计算 PSD 并将其存储在 pxx 的对应列中。如果 x 是实数值，则 pxx 是单边 PSD 估计值。如果 x 是复数值，则 pxx 是双边 PSD 估计值。默认情况下，x 被分成尽可能长的段，以获得接近但不超出 8 个重叠率为 50% 的段。对每段使用一个汉明窗进行加窗处理。对各个修正周期图求平均值以获得 PSD 估计值。如果您无法将 x 的长度精确地分成重叠率为 50% 的整数段数，x 会相应地截断。
[___, 在 pxxc] = pwelch(___,'ConfidenceLevel',probability)pxxc 中返回 PSD 估计值的 probability × 100% 置信区间。
[___] = pwelch(___, 在 spectrumtype)spectrumtype 指定为 'psd' 时返回 PSD 估计值，在 spectrumtype 指定为 'power' 时返回功率谱。
不带输出参量的 pwelch(___) 在当前图窗窗口中绘制韦尔奇 PSD 估计值。
示例
获取一个输入信号的韦尔奇 PSD 估计值，该输入信号由角频率为 弧度/采样点的离散时间正弦信号和加性 白噪声组成。
创建一个角频率为 弧度/采样点并具有加性 白噪声的正弦波。重置随机数生成器以获得可重现的结果。信号的长度为 个采样。
rng default
n = 0:319;
x = cos(pi/4*n)+randn(size(n));使用默认汉明窗和 DFT 长度获取韦尔奇 PSD 估计值。默认段长度为 71 个采样，DFT 长度为 256 个点，产生 弧度/采样点的频率分辨率。由于信号是实数值信号，周期图为单边型并且有 256/2+1 个点。绘制韦尔奇 PSD 估计值。
pxx = pwelch(x); pwelch(x)
重复该计算。
将信号分成长度为 的若干段。此操作等效于将信号分成尽可能长的段，以获得接近但不超出 8 个叠加率为 50% 的段。
使用汉明窗对这些段进行加窗处理。
指定连续段之间的重叠率为 50%
为了计算 FFT，请使用 个点，其中 。
验证这两种方法是否给出相同的结果。
Nx = length(x); nsc = floor(Nx/4.5); nov = floor(nsc/2); nff = max(256,2^nextpow2(nsc)); t = pwelch(x,hamming(nsc),nov,nff); maxerr = max(abs(abs(t(:))-abs(pxx(:))))
maxerr = 0
将信号分成 8 个等长的段，段之间的重叠率为 50%。指定与上一步相同的 FFT 长度。计算韦尔奇 PSD 估计值，并验证它是否与前两个过程的结果相同。
ns = 8; ov = 0.5; lsc = floor(Nx/(ns-(ns-1)*ov)); t = pwelch(x,lsc,floor(ov*lsc),nff); maxerr = max(abs(abs(t(:))-abs(pxx(:))))
maxerr = 0
获取一个输入信号的韦尔奇 PSD 估计值，该输入信号由角频率为 弧度/采样点的离散时间正弦信号和加性 白噪声组成。
创建一个角频率为 弧度/采样点并具有加性 白噪声的正弦波。重置随机数生成器以获得可重现的结果。信号有 512 个采样。
rng default
n = 0:511;
x = cos(pi/3*n)+randn(size(n));通过将信号分成长度为 132 个采样的段来获取韦尔奇 PSD 估计值。信号段乘以长度为 132 个采样的汉明窗。未指定重叠采样数，因此设置为 132/2 = 66。DFT 长度为 256 个点，产生 弧度/采样点的频率分辨率。由于信号是实数值信号，PSD 估计值为单边型并且有 256/2+1 = 129 个点。将 PSD 绘制为归一化频率的函数。
segmentLength = 132;
[pxx,w] = pwelch(x,segmentLength);
plot(w/pi,10*log10(pxx))
xlabel('\omega / \pi')
获取一个输入信号的韦尔奇 PSD 估计值，该输入信号由角频率为 弧度/采样点的离散时间正弦信号和加性 白噪声组成。
创建一个角频率为 弧度/采样点并具有加性 白噪声的正弦波。重置随机数生成器以获得可重现的结果。信号长度为 320 个采样。
rng default
n = 0:319;
x = cos(pi/4*n)+randn(size(n));通过将信号分成长度为 100 个采样的段来获取韦尔奇 PSD 估计值。信号段乘以长度为 100 个采样的汉明窗。重叠采样数为 25。DFT 长度为 256 个点，产生 弧度/采样点的频率分辨率。由于信号是实数值信号，PSD 估计值为单边型并且有 256/2+1 个点。
segmentLength = 100; noverlap = 25; pxx = pwelch(x,segmentLength,noverlap); plot(10*log10(pxx))
获取一个输入信号的韦尔奇 PSD 估计值，该输入信号由角频率为 弧度/采样点的离散时间正弦信号和加性 白噪声组成。
创建一个角频率为 弧度/采样点并具有加性 白噪声的正弦波。重置随机数生成器以获得可重现的结果。信号长度为 320 个采样。
rng default
n = 0:319;
x = cos(pi/4*n) + randn(size(n));通过将信号分成长度为 100 个采样的段来获取韦尔奇 PSD 估计值。使用默认重叠率 50%。指定 DFT 长度为 640 个点，以使 弧度/采样点的频率对应于一个 DFT bin (bin 81)。由于信号是实数值信号，PSD 估计值为单边型并且有 640/2+1 个点。
segmentLength = 100; nfft = 640; pxx = pwelch(x,segmentLength,[],nfft); plot(10*log10(pxx)) xlabel('rad/sample') ylabel('dB / (rad/sample)')
创建一个信号，该信号包含在加性 N(0,1) 白噪声中的一个 100 Hz 正弦波。重置随机数生成器以获得可重现的结果。采样率为 1 kHz，信号持续时间为 5 秒。
rng default
fs = 1000;
t = 0:1/fs:5-1/fs;
x = cos(2*pi*100*t) + randn(size(t));获取上述信号的韦尔奇交叠分段平均法 PSD 估计值。使用一个长度为 500 个采样的段，其中有 300 个重叠采样。使用 500 个 DFT 点，以使 100 Hz 恰好位于一个 DFT bin 上。输入采样率以输出以 Hz 为单位的频率向量。绘制结果。
[pxx,f] = pwelch(x,500,300,500,fs); plot(f,10*log10(pxx)) xlabel('Frequency (Hz)') ylabel('PSD (dB/Hz)')
创建一个信号，该信号由三个含噪正弦波和一个以 200 kHz 进行 0.1 秒采样的啁啾组成。这些正弦波的频率分别为 1 kHz、10 kHz 和 20 kHz。这些正弦波具有不同振幅和噪声电平。无噪声啁啾的频率从 20 kHz 开始，并在采样期间线性增大到 30 kHz。
Fs = 200e3; Fc = [1 10 20]'*1e3; Ns = 0.1*Fs; t = (0:Ns-1)/Fs; x = [1 1/10 10]*sin(2*pi*Fc*t)+[1/200 1/2000 1/20]*randn(3,Ns); x = x+chirp(t,20e3,t(end),30e3);
计算信号的韦尔奇 PSD 估计值以及最大值保持和最小值保持频谱。绘制结果。
[pxx,f] = pwelch(x,[],[],[],Fs); pmax = pwelch(x,[],[],[],Fs,'maxhold'); pmin = pwelch(x,[],[],[],Fs,'minhold'); plot(f,pow2db(pxx)) hold on plot(f,pow2db([pmax pmin]),':') hold off xlabel('Frequency (Hz)') ylabel('PSD (dB/Hz)') legend('pwelch','maxhold','minhold')
重复此过程，这次计算中心化功率谱估计值。
[pxx,f] = pwelch(x,[],[],[],Fs,'centered','power'); pmax = pwelch(x,[],[],[],Fs,'maxhold','centered','power'); pmin = pwelch(x,[],[],[],Fs,'minhold','centered','power'); plot(f,pow2db(pxx)) hold on plot(f,pow2db([pmax pmin]),':') hold off xlabel('Frequency (Hz)') ylabel('Power (dB)') legend('pwelch','maxhold','minhold')
此示例说明置信边界与韦尔奇交叠分段平均法 (WOSA) PSD 估计值的用法。在韦尔奇估计值中，如果某个频率的置信边界下限超出周围 PSD 估计值的置信边界上限，则清楚地表明时间序列中存在显著的振荡，虽然这不是统计意义的必要条件。
创建一个信号，该信号包含在加性白 N(0,1) 噪声中叠加的 100 Hz 和 150 Hz 正弦波。这两个正弦波的振幅为 1。采样率为 1 kHz。重置随机数生成器以获得可重现的结果。
rng default
fs = 1000;
t = 0:1/fs:1-1/fs;
x = cos(2*pi*100*t)+sin(2*pi*150*t)+randn(size(t));获取具有 95% 置信边界的 WOSA 估计值。将段长度设置为 200 个采样，并使段重叠 50%（100 个采样）。绘制 WOSA PSD 估计值以及置信区间，并放大在 100 Hz 和 150 Hz 附近的感兴趣频率区域。
L = 200; noverlap = 100; [pxx,f,pxxc] = pwelch(x,hamming(L),noverlap,200,fs,... 'ConfidenceLevel',0.95); plot(f,10*log10(pxx)) hold on plot(f,10*log10(pxxc),'-.') hold off xlim([25 250]) xlabel('Frequency (Hz)') ylabel('PSD (dB/Hz)') title('Welch Estimate with 95%-Confidence Bounds')
紧靠 100 Hz 和 150 Hz 的置信边界下限显著高于 100 Hz 和 150 Hz 以外的置信边界上限。
创建一个信号，该信号包含在加性 白噪声中的一个 100 Hz 正弦波。重置随机数生成器以获得可重现的结果。采样率为 1 kHz，信号持续时间为 5 秒。
rng default
fs = 1000;
t = 0:1/fs:5-1/fs;
noisevar = 1/4;
x = cos(2*pi*100*t)+sqrt(noisevar)*randn(size(t));使用韦尔奇方法获取以 DC 为中心的功率谱。使用 500 个采样（具有 300 个重叠采样）的段长度和 500 个点的 DFT 长度。绘制结果。
[pxx,f] = pwelch(x,500,300,500,fs,'centered','power'); plot(f,10*log10(pxx)) xlabel('Frequency (Hz)') ylabel('Magnitude (dB)') grid
您可以看到在 -100 Hz 和 100 Hz 处的功率接近于 1/4 的预期功率，这是振幅为 1 的实数值正弦波的功率。与 1/4 的偏差是加性噪声造成的。
生成一个多通道信号的 1024 个采样，该信号包含在加性 高斯白噪声中叠加的三个正弦波。这些正弦波的频率为 、 和 弧度/采样点。使用韦尔奇方法估计信号的 PSD 并对其绘图。
N = 1024; n = 0:N-1; w = pi./[2;3;4]; x = cos(w*n)' + randn(length(n),3); pwelch(x)
输入参数
输出参量
详细信息
参考
[1] Hayes, Monson H. Statistical Digital Signal Processing and Modeling. New York: John Wiley & Sons, 1996.
[2] Stoica, Petre, and Randolph Moses. Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice Hall, 2005.
