matlab - Mathcad 到 Matlab - 白噪声、fft 和 NPS 测试

Question

我正在尝试在 Matlab 中编写简单的函数来计算和绘制噪声功率谱 (NPS)。首先，我想测试一下我从老师那里得到的算法是否正常。这是它(它是用 mathcad 制作的)

所以我尝试将其复制粘贴到 Matlab 脚本中，最终得到以下代码:

clear all;
clc;


N=1000;
O=1024;
mn=zeros(N,O);
n0=1500;
s=sqrt(n0);
W=zeros(N,O/2);
W1=zeros(N,O);

for k=1:N
    for l=1:O
        mn(k,l)=n0+round(sin(randn)*s);
    end
end

for k=1:N
    for l=1:O
        mn(k,l)=mn(k,l)-n0;
    end
end

for k=1:N
    W1(k,:)=fft(mn(k,:));
end

for k=1:N
   for l=1:O/2
       W(k,l)=W1(k,l);
   end
end

NPS1=(abs(W)).^2;
NPS2=sum(NPS1);
NPS=(1/N)*NPS2;

plot(NPS);

我没有使用泊松分布，并且我已经切换了行列索引，但这应该不重要(对吗？)。问题是我图中的值几乎比预期大 400 倍。

它应该是这样的:

我试图找出我做错了什么，但经过相当长一段时间的测试一些更改后，我又回到了原点...唯一让我担心的是，也许 Matlab fft 函数的工作方式与Mathcad 中使用的一个(不能说我完全理解它)。 哪位好心人能告诉我这是否与 fft 函数有关？或者我只是一个盲目的傻瓜，看不到他所犯的愚蠢错误？ 干杯，抱歉我的英语不好。

[编辑]

过了一段时间后，我的老师要求我检查此方法是否适用于相关(有点)噪声，因为它(再次)适用于 mathcad。 关联后，其 NPS 应以更高的频率“下降”。问题是事实并非如此。我用来测试的代码如下所示:

clear all;
clc;

N=1000;

mn = poissrnd(N, N, N);
dataw=zeros(N);

for k=1:N ## loop used for my teacher's correlation method
    for l=1:N
        if l>1 && l<N
            dataw(k,l)=dataw(k,l)+mn(k,l)*0.5+mn(k,l-1)*0.25+mn(k,l+1)*0.25;
        elseif l==1
            dataw(k,l)=dataw(k,l)+mn(k,l)*0.75+mn(k,l+1)*0.25;
        else
            dataw(k,l)=dataw(k,l)+mn(k,l)*0.75+mn(k,l-1)*0.25;
        end
    end
end

dataw = dataw - mean(dataw(:));
W1 = (1/sqrt(N))*fft(dataw, [], 1);

NPS1=(abs(W1)).^2;
NPS2=sum(NPS1);
NPS=(1/N)*NPS2;

plot(NPS);

我对 rayryeng 修复的代码所做的唯一更改是将噪声矩阵设置为平方 (1000x1000)，将平均值设置为 1000，并使用整个变换后的向量 W1 而不是其“一半”。 我知道它对我的老师有用，但对我不起作用...关于 matlab fft 还有其他我忽略的东西还是我使用的“相关方法”？

经过几次快速更改后，在 Mathcad 中添加了它的外观(有一些细微的差异，但总的来说，它显示了我应该获得的效果)。它切断了扫描的开头，但它与我在本文开头放置的 1 完全相同。

[编辑2]

Nvm，这只是fft函数中的维度问题。改成fft(dataw, [], 2)后看起来更好了。

Answer 1

它不起作用的主要原因是MathCad和MATLAB之间FFT的比例因子。对于 MathCad，有一个额外比例因子 1/sqrt(N)，而 MATLAB不包含该比例因子。因此，如果您想模拟使用 MathCad 看到的结果，则需要将 FFT 结果乘以该比例因子。

另外，我对您的代码有一些建议:

1. 我们可以将其完全矢量化，无需任何循环
2. fft 和 randn 等函数可以对矩阵进行运算，您可以将该函数专门应用于某一特定维度。

请注意，我已将您的随机噪声分布替换为泊松随机噪声(来自 poissrnd )，以便我可以模仿您老师看到的结果。

---

本质上，您的代码可以替换为:

clear all;
clc;

N=1000;
O=1024;
n0=1500;
s=sqrt(n0);

%mn = round(sin(randn(N,O)*s));
mn = poissrnd(n0, N, O); %// CHANGE
mn = mn - mean(mn(:)); %// Remove mean
W1 = (1/sqrt(N))*fft(mn, [], 1); %// CHANGE FROM ABOVE
W = W1(:,1:O/2);

NPS1=(abs(W)).^2;
NPS2=sum(NPS1);
NPS=(1/N)*NPS2;

plot(NPS);

---

请注意，您在生成随机数据时添加了平均值 1500....只是再次从中减去 1500，而不对偏移数据进行任何处理。我刚刚从正弦舍入随机噪声的代码中删除了它。我已将该代码注释掉，因为无论如何我现在都不会运行该代码。另请注意，randn 可以获取行数和列数，以便您可以生成随机值矩阵。此外，fft 可以对行或列进行运算，并将该维度中的每个信号视为一维信号。在本例中，您希望对每一列进行操作，对行进行处理，这就是为什么我们将1参数指定为第三个参数。

这是我运行上面的代码时得到的结果:

---

您会看到它徘徊在 1500 均值左右，这正是我们从具有 lambda=1500 的泊松随机分布中得出的结果。如果您确实想让图表看起来像您老师的图表，请将 y 轴的限制从 0 更改为 2000，如下所示:

ylim([0 2000]);

因此我们得到: