Matlab - 创建具有正确缩放比例的 psd

Question

我想分析一个音频数据（.wav，pcm，32k 作为采样率）并使用轴 Sxx（瓦特/赫兹而不是 db）和 f（赫兹）创建它的 psd。

所以我会先读出音频数据：

[x,fs]=audioread('test.wav');

在这之后我遇到了一些问题，因为我真的不知道如何进行，而且 Matlab 总是告诉我将来不支持 psd 函数，我应该使用 pwelch..（也尝试构建autocorr，然后使用傅立叶到达 Sxx，但效果并不好）

那么谁能告诉我如何从我的向量 x 得到一个以瓦特/赫兹为单位的 psd 值的向量，然后再绘制它？

非常感谢您的各种帮助！ :)

Update1：是的，我确实阅读了 pwelch 的文档，但恐怕我的英语太糟糕了，无法完全理解。 所以如果我使用 psd 文档：

nfft = 2^nextpow2(length(x));
Pxx = abs(fft(x,nfft)).^2/length(x)/fs;
Hpsd = dspdata.psd(Pxx(1:length(Pxx)/2),'fs',fs);  
plot(Hpsd) 

我能够以正确的频率获得峰值在 db 中的图。 （虽然我不知道 dspdata.psd 是如何工作的）

我试过了：

[Pyy,f]=pwelch(x,fs)
plot(Pyy)

这给了我一个非分贝规模，但峰值的频率错误

更新 2： 首先非常感谢您的详细解答！目前我正在学习我的 matlabskills 以及我的英语，但所有具体的技术术语都让我很难过.. 当在明确频率为 1khz 的 wav 数据上使用 pwelch 示例时，该图向我显示了大约 0.14 附近的峰值，它可能仍然是一个特殊比例的 x 轴吗？

如果我这样尝试：

[y,fs]=audioread('test.wav');
N=length(y);
bin_vals=0:N-1;
fax_Hz= bin_vals*fs/N;
N_2=ceil(N/2);
Y=fft(y);
pyy=Y.*conj(Y);
plot(fax_Hz(1:N_2),pyy(1:N_2))

结果似乎正确（这种方式正确吗？），但我仍然需要一些时间来寻找以 W/Hz 为单位显示 y 轴的正确方法，因为我不知道音频信号是如何创建的。

更新 3： http://s000.tinyupload.com/index.php?file_id=33803229773204653857 此 wav 文件应具有 1khz 的主频率，持续时间为 3 秒，采样频率为 44100Hz。 （如果我绘制从 audioread 接收到的数据，振荡似乎是合理的）

与

[y,fs]=audioread('1khz.wav');
[pyy,f]=pwelch(y,fs);
plot(f,pyy)

我在 x 轴上的 0.14 处得到一个峰值。

如果我使用

[y,fs]=audioread('1khz.wav');
[pyy,f]=pwelch(y,[],[],[],fs);
plot(f,pyy)

相反，峰值在 1000 处。这样对吗？我如何解释 y 轴上的差异缩放？ （pwelch 与 abs 的平方） 我还想问是否有可能在matlab中获得awgn的平面psd？ （因为你只有有限的元素，我不知道如何到达那里）

再次感谢您的详细支持！

更新 4 @A.唐达 所以我有一个新问题，我认为可能有必要更详细地讨论它。所以我的计划基本上是做到以下几点：

1. 读取音频数据 ([y,fs]) 并生成具有一定 SNR ([n,fs]) 的白噪声
2. 生成与 PSD(n) 相似的 PSD(y) 形状的滤波器 H
3. 生成一个反向过滤器 G=H^(-1) 来恢复 H 的效果。

我的问题是，使用 pwelch，得到的 pyy 的向量长度远小于 y 的向量长度。由于我的过滤器由 P=sqrt(pnn/pyy) 确定，因此我无法乘以 fft(y)*H，因此没有得到任何结果。 你知道这个问题有什么帮助吗？ 或者有没有办法从 PSD（Welch 估计）返回到正常信号（如 pwelch 的反函数）？

Answer 1

在 psd 文档中的示例中，您自己计算 psd 估计值，然后将其放入 dspdata.psd 容器并绘制它。 dspdata.psd 数据在这里为您所做的基本上是计算频率轴并将其提供给 plot 命令，仅此而已。您会得到频谱密度估计图，但这是您使用fft 自己计算的图，这是您可以获得的最简单和最差的 psd 估计值，即所谓的周期图。

您对pwelch 的使用几乎是正确的，您只是忘记在绘图中使用频率轴信息。

[Pyy,f]=pwelch(x,fs)
plot(f,Pyy)

应该以正确的频率为您提供峰值。

您对pwelch 的使用几乎是正确的，但您必须将采样频率作为第 5 个参数，然后在绘图中使用频率轴信息。

[Pyy,f]=pwelch(y,[],[],[],fs);
plot(f,Pyy)

应该给你正确频率的峰值。

pwelch 为您提供的是信号在 Hz 上的频谱密度。因此，正确的轴标签将是

xlabel('frequency (Hz)')
ylabel('psd (1/Hz)')

您给pwelch 的信号是一个没有物理尺寸的纯数字序列。通过指定采样率，时间轴得到一个物理单位，s，因此得到的频率以 Hz 为单位，密度以 1/Hz 为单位。但是您的时间序列值仍然没有物理维度，因此密度不能与 W 之类的东西相关。您的音频信号是通过校准的 A/D 转换器获得的吗？如果是，您应该能够将您的数据与物理维度和单位相关联，但这是一个重要的步骤。

就个人而言，我真的建议你复习一下你的英语，因为在没有正确理解文档的情况下使用软件，尤其是编程接口，会导致灾难。