Matlab - 创建具有正确缩放比例的 psd
Matlab - creating psd with right scaling
我想分析音频数据（.wav 与 pcm，32k 作为采样率）并使用轴 Sxx（watts/hertz 不是 db）和 f（赫兹）创建它的 psd。
所以我会先读出音频数据：
[x,fs]=audioread('test.wav');
在此之后我遇到了一些问题，因为我真的不知道如何继续，而且 Matlab 总是告诉我未来将不支持 psd 函数，我应该使用 pwelch..（也尝试过构建 autocorr 然后使用 fourier 得到 Sxx 但效果不是很好）
那么谁能告诉我如何从我的向量 x 得到一个具有 watts/hertz 中的 psdvalues 的向量，然后绘制它？
非常感谢各种帮助！ :)
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)
我能够以正确的频率获得峰值的分贝图。 （虽然我不知道 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.Donda
所以我有一个新问题，我认为可能有必要更详细地讨论一下。所以我的计划基本上是做以下事情：
- 读取音频数据 ([y,fs]) 并生成具有特定 SNR ([n,fs]) 的白噪声
- 生成滤波器 H，使 PSD(y) 的形状类似于 PSD(n)
- 生成一个反向过滤器 G=H^(-1) 以恢复 H 的效果。
我的问题是使用 pwelch 时，pyy 的矢量长度比 y 的矢量长度小得多。由于我的过滤器由 P=sqrt(pnn/pyy) 决定，我不能乘以 fft(y)*H，因此得不到任何结果。
你知道这个问题有什么帮助吗？
或者有没有办法从 PSD（Welch 估计）返回到正常信号（如 pwelch 的反函数）？
在 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 转换器获得？如果是，您应该能够将您的数据与物理尺寸和单位相关联，但这是一个重要的步骤。
就个人而言，我真的建议您复习一下英语，因为在没有正确理解文档的情况下使用软件，尤其是编程界面，会导致灾难。
我想分析音频数据（.wav 与 pcm，32k 作为采样率）并使用轴 Sxx（watts/hertz 不是 db）和 f（赫兹）创建它的 psd。
所以我会先读出音频数据：
[x,fs]=audioread('test.wav');
在此之后我遇到了一些问题，因为我真的不知道如何继续，而且 Matlab 总是告诉我未来将不支持 psd 函数，我应该使用 pwelch..（也尝试过构建 autocorr 然后使用 fourier 得到 Sxx 但效果不是很好）
那么谁能告诉我如何从我的向量 x 得到一个具有 watts/hertz 中的 psdvalues 的向量，然后绘制它？
非常感谢各种帮助！ :)
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)
我能够以正确的频率获得峰值的分贝图。 （虽然我不知道 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.Donda 所以我有一个新问题，我认为可能有必要更详细地讨论一下。所以我的计划基本上是做以下事情：
- 读取音频数据 ([y,fs]) 并生成具有特定 SNR ([n,fs]) 的白噪声
- 生成滤波器 H，使 PSD(y) 的形状类似于 PSD(n)
- 生成一个反向过滤器 G=H^(-1) 以恢复 H 的效果。
我的问题是使用 pwelch 时，pyy 的矢量长度比 y 的矢量长度小得多。由于我的过滤器由 P=sqrt(pnn/pyy) 决定，我不能乘以 fft(y)*H，因此得不到任何结果。 你知道这个问题有什么帮助吗？ 或者有没有办法从 PSD（Welch 估计）返回到正常信号（如 pwelch 的反函数）？
在 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 转换器获得？如果是，您应该能够将您的数据与物理尺寸和单位相关联，但这是一个重要的步骤。
就个人而言，我真的建议您复习一下英语，因为在没有正确理解文档的情况下使用软件，尤其是编程界面，会导致灾难。