2

我正在尝试在 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) 后看起来更好。

4

1 回答 1

2

它不起作用的主要原因是由于 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]);

因此我们得到：

在此处输入图像描述

于 2014-12-12T18:30:42.553 回答