0

我正在尝试在 matlab 中测量随机过程的 PSD，但我不知道该怎么做。我在这里发布了完全相同的问题，但我想我可能在这里有更多的运气。

随机过程描述风速，并由实数向量表示。每个条目对应于空间中某个点的风速，以 m/s 为单位测量。这些点相距 0.0005 m。如何测量和绘制 PSD？我们称之为向量V。我的第一个想法是使用

[p, w] = pwelch(V);
loglog(w,p);

但这是正确的吗？问题是，我得到了一个解析表达式，PSD 应该（理论上）匹配。通过用这两行代码绘制它，它看起来全错了。具体来说，它看起来好像需要平移和缩放。除此之外，两者的形状相似。

在此处输入图像描述

更新： 上面的图像实际上并没有描述通过使用pwelch单个向量获得的 PSD，而是 200 个向量的 PSD 的平均值，因为这些向量源于数值模拟。正如建议的那样，我尝试按 2*pi/0.0005 进行缩放。我看到您实际上可以将此信息提供给pwelch. 所以我尝试使用代码

[p, w] = pwelch(V,[],[],[],2*pi/0.0005);
loglog(w,p);

反而。如下所示，它看起来好多了。然而，它仍然不完美。这只是我应该期待的吗？顺便说一句，取平方根不是答案。不过谢谢你的建议。一方面，它应该遵循 Kolmogorov 的 -5/3 定律，它现在确实如此（它遵循绿线，斜率为 -5/3）。我试图与之匹配的函数是 Shkarofsky 谱密度函数，它是 Shkarofsky 相关函数的一维傅里叶变换。不能在网站上标记数学吗？

在此处输入图像描述

更新2： 我已经尝试[p, w] = pwelch(V,[],[],[],1/0.0005);按照我的建议使用。但正如你所见，它仍然不完全匹配。我很难准确地解释我在寻找什么。但我想要（或者，我所期望的）是计算和分析 PSD 的下降同时发生，并以相同的速度下降。数据来自湍流模拟。解析表达式已拟合到湍流的实际测量值，其中也存在这种下降。我根本不是专家，但据我所知，下降发生在较小的长度尺度上，因为由于空气的粘性，能量会消散。

在此处输入图像描述

4

3 回答 3

1

使用标准方程来获得 aPSD怎么样？我会这样做：

Sxx(f) = (fft(x(t)).*conj(fft(x(t))))*(dt^2);

或者

Sxx = fftshift(abs(fft(x(t))))*(dt^2);

然后，如果您真的需要，您可能会考虑应用一个窗口标准，例如

  • 汉宁
  • 汉明
  • 韦尔奇

这只会以某种方式过滤您的PSD.

于 2013-08-15T18:35:32.520 回答
0

大概您需要将频率（波数）重新调整为 1/m 的单位。pwelch应该重新调整频率单位，因为正如文档解释的那样：

W 是估计 PSD 的归一化频率向量。W 的单位为 rad/sample。

即兴发挥我的猜测是比例因子是

 scale = 1/0.0005/(2*pi);

或 318.3 (m^-1)。

至于强度，看起来取平方根可能会有所帮助。也许您的方程报告的是强度，而不是 PSD？

编辑

正如您所指出的，由于分析和计算的 PSD 具有几乎相同的斜率，因此它们似乎遵循高达 800 m^-1 的相似幂律。我不确定您需要指数或偏移量到什么程度才能满足特定模型，而且我不熟悉这个特定理论。

至于高波数的明显不一致，我要指出您正在进入非常小的数字域，因此（1）浮点问题和（2）噪声可能潜伏着。另一方面，在计算的 PSD 中看起来非常漂亮的下降看起来非常真实，但我对此没有任何解释（也许你的噪音不是白色的？）。

您可能想在 matlab 中心查看此提交，因为它可能有用。

编辑#2

在检查 的文档后pwelch，您似乎应该通过1/0.0005（采样率）而不是2*pi/0.0005. 这不应该影响斜率，但会影响截距。

于 2013-08-15T15:28:48.033 回答
0

模拟结果中 PSD 的下降看起来类似于我在使用低阶方法对原始数据进行插值时在数据中看到的混叠伪影。为了让这一点更清楚 - 假设我的原始数据间隔为 0.002m，但在清理丢失数据的过程中，试图节省空间，无论如何，我将这些数据线性插值到 0.005m 间隔。线性插值的频率响应表现不佳，并且会在频谱的高波数端引入波峰和波谷。

光谱估计有不同的约定——波数单位是 1/m 还是弧度/m。单面光谱或双面光谱。

帮助

显示默认设置返回单边频谱，即某个频率 ω 的 bin 将包括 +ω 和 -ω 的功率密度。您应该仔细检查您要比较的理想化光谱是否也是单面光谱。否则，您需要将一侧光谱的值减半，以获得代表两侧光谱的 +ω 侧的值。

我同意 Try Hard 的观点，即应该指定为 pwelch 的循环频率（通常为 Hz，或在本例中为 1/m）。也就是说，从 pwelch 返回的频率向量也在这些单位中。分析光谱公式通常以角频率形式编写，因此您需要确保以弧度/m 的形式对其进行评估，但要按比例缩小到 1/m 进行绘图。

于 2014-07-11T14:46:06.253 回答