MATLAB中功率谱密度计算pwelch函数使用详解_C/C++

window--------窗口，指定为行向量或列向量或整数。如果 window 是一个向量，pwelch 将 x 划分为长度等于 window 长度的重叠段，然后将每个信号段与 window 中指定的向量相乘。如果window是整数，则将pwelch分成长度等于整数值的段，并使用等长的汉明窗。如果x的长度不能精确地划分为具有noverlap数量的重叠样本的整数段，则x被相应地截断。如果将 window 指定为空，则使用默认的 hamming 窗来获取 x 的 8 段，其中具有 noverlap 重叠样本。

noverlap----------重叠样本的数量，指定为小于窗口长度的正整数。如果省略 noverlap 或 noverlap 指定为空，则使用一个值来获得段之间 50% 的重叠，即默认50%重叠。

nfft--------dft 点数，指定为正整数。对于实值输入信号 x（psd 估计值），如果 nfft 为偶数，则 pxx 的长度为 (nfft/2 + 1)；如果 nfft 为奇数，则 pxx 的长度为 (nfft + 1)/2。对于复值输入信号 x，psd 估计的长度始终为 nfft。如果 nfft 指定为空，则使用默认的 nfft。如果 nfft 大于window长度，则数据用零填充。如果 nfft 小于window长度，则使用 datawrap包装该段，使长度等于nfft。建议两者相等。

fs-------采样频率，指定为正标量。采样率是单位时间内的采样数。如果时间单位是秒，那么采样率的单位是hz。

freqrange--------------psd 估计的频率范围，指定为“单边”、“双边”或“中心”之一。对于实值信号，默认值为“单侧”；对于复值信号，默认值为“双侧”。每个选项对应的频率范围为

'oneside' — 返回实值输入信号 x 的单侧 psd 估计。如果 nfft 为偶数，则 pxx 的长度为 nfft/2 + 1，如果 nfft 为奇数，则 pxx 的长度为 (nfft + 1)/2，间隔为 [0,π) rad/sample。当 fs 可选指定时，对于偶数和奇数长度 nfft，相应的间隔分别为 [0,fs/2] 周期/单位时间和 [0,fs/2) 周期/单位时间。该函数将除 0 和奈奎斯特频率之外的所有频率的功率乘以 2，以保持总功率不变。

'twoside' - 返回实值或复值输入 x 的两侧 psd 估计值。在这种情况下，pxx 的长度为 nfft，并在区间 [0,2π) rad/sample 上计算。当指定fs时，间隔为[0,fs)周期/单位时间。

'centered' - 返回实值或复值输入 x 的中心两侧 psd 估计值。在这种情况下，pxx 的长度为 nfft，并且在偶数长度 nfft 的间隔 (–π,π] rad/sample 和奇数长度 nfft 的 (–π,π) rad/sample 区间内计算。当 fs 可选地指定时，相应的对于偶数和奇数长度 nfft，间隔分别为 (–fs/2, fs/2] 周期/单位时间和 (–fs/2, fs/2) 周期/单位时间。

spectrumtype------------功率谱类型，指定为“psd”或“power”之一。默认“psd”，将返回功率谱密度。指定“power”可通过窗口的等效噪声带宽来缩放 psd 的每个估计值。使用“power”选项可以获得每个频率的功率估计值。

三、pxx = pwelch(x)示例

创建一个角频率为π/4 rad的正弦波，外加n(0,1)白噪声。信号的长度是320个样本。得到其welch psd估计。

clc;
clear;
close all;

rng default

n = 0:319;
x = cos(pi/4*n)+randn(size(n));

pxx = pwelch(x);

可见，默认横轴是归一化的角频率，纵轴是取了10log10( )的db功率谱密度。

关于归一化频率，参考：滤波器设计中的频率归一化问题_滤波器归一化频率-csdn博客

解释如下：

信号处理工具箱中经常使用的频率是nyquist频率，它被定义为采样频率的一半，在滤波器的结束选择和设计当中的截止频率均使用nyquist频率进行归一化处理。

例如，对于一个采样频率为1000hz的系统，300hz的归一化即为300/500=0.6。归一化频率的范围在[0,1]之间。如果要将归一化频率转换为角频率，则将归一化频率乘以pi；如果将归一化频率转换成hz，则将归一化频率乘以采样频率的一半。

采样率的一半是最高频率，认为是１，那么真实频率和最高频率的比值就是归一化频率！也叫数字频率。
将信号分成长度为nsc=⌊nx/4.5⌋。这个动作相当于将信号分成尽可能长的段，以获得接近但不超过8个重叠50%的段。使用汉明窗口显示各部分。指定相邻部分之间50%的重叠

要计算fft，使用max(256,2^p)，其中p=[log2nsc⌉。

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(:))))

默认分成8段，每段的长度为nx/4.5=320/4.5=71.1111，舍去多余的数据，分段长度为71，并用同等长度的汉明窗；重叠长度为50%，则nov=floor(71/2)=35；dft的点数取每段长度最接近的2的整数次幂和256的最大值，最接近的2的整数次幂是比每段长度长的最接近的2的整数次幂，所以dft计算的时候如果dft点数大于每段长度。会自动补0。

对于实值信号，默认值为“单侧”psd，所以计算dft点数为256，估计的psd长度只有256/2+1=129点的长度。

四、[pxx,f]=pwelch(x,window,noverlap,nfft,fs)示例

创建一个由100hz正弦波和加性n(0,1)白噪声组成的信号。采样率为1khz，信号持续时间为5秒。使用500个样本和300个重叠样本的段长度，使用500个dft点。

clc;
clear;
close all;

rng default

fs = 1000;
t = 0:1/fs:5-1/fs;
x = cos(2*pi*100*t) + randn(size(t));

[pxx,f] = pwelch(x,500,300,500,fs);

plot(f,10*log10(pxx))

xlabel('frequency (hz)')
ylabel('psd (db/hz)')

sum(abs(x).^2)/length(x) %信号功率
p = fs/2/length(pxx) * (sum(pxx) - 0.5*(pxx(1) + pxx(end))) %梯形法积分

验证功率，对功率谱密度积分发现和信号功率几乎相等。

增大点数

clc;
clear;
close all;

rng default

fs = 10000;
t = 0:1/fs:5-1/fs;
x = cos(2*pi*100*t) + randn(size(t));

[pxx,f] = pwelch(x,5000,3000,5000,fs);

plot(f,10*log10(pxx))

xlabel('frequency (hz)')
ylabel('psd (db/hz)')

sum(abs(x).^2)/length(x) %信号功率
p = fs/2/length(pxx) * (sum(pxx) - 0.5*(pxx(1) + pxx(end))) %梯形法积分

增大段长度，会发现功率谱估计的准，功率谱密度积分的值更接近信号的功率。

四、[pxx,f] = pwelch(x,window,noverlap,nfft,fs,freqrange,spectrumtype)示例

创建一个由100hz正弦波和加性n(0,1)白噪声组成的信号。采样率为1khz，信号持续时间为5秒。使用500个样本和300个重叠样本的段长度，使用500个dft点。

clc;
clear;
close all;

rng default

fs = 1000;
t = 0:1/fs:10-1/fs;

x = cos(2*pi*100*t)+randn(size(t))+1;

[pxx,f] = pwelch(x,500,300,500,fs,'centered');

plot(f,10*log10(pxx))
xlabel('frequency (hz)')
ylabel('psd (db/hz)')
e=sum(abs(x).^2)  %信号能量
p=sum(abs(x).^2)/length(x) %信号功率
p = fs/length(pxx) * (sum(pxx) - 0.5*(pxx(1) + pxx(end))) %梯形法积分�

[pxx,f] = pwelch(x,500,300,500,fs,'centered','power');
figure(2)
plot(f,10*log10(pxx))
xlabel('frequency (hz)')
ylabel('magnitude (db)')
grid

使用参数centered得到双边功率谱密度：

使用参数power得到双边功率谱（不是双边功率谱密度）：

观察图可知，与信号设置的功率吻合。

取对数画出如下：

可见，将功率谱密度积分和信号功率相等；观察功率谱图可知，每个频率对应的点显示了该频点的功率大小。

五、多通道功率谱估计

在加性n(0,1)高斯白噪声中生成由三个正弦波组成的多通道信号的1024个样本。正弦波的频率分别为100、200、300hz，采样频率为1000hz。使用welch的方法估计信号的psd并绘制它。

clc;
clear;
close all;

rng default

fs = 1000;
t = 0:1/fs:5-1/fs;
f = [100;200;300];
x = cos(2*pi*f*t)' + randn(length(t),3);

[pxx,f] = pwelch(x,500,300,500,fs);

plot(f,10*log10(pxx));

xlabel('frequency (hz)')
ylabel('psd (db/hz)')