Matlab中的Mathcad白噪声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);

我不使用Poisson分布，我已经切换了行列索引，但这不重要(对吗？)问题是，我的作品中的值几乎比预期的要大400倍。

下面是它应该是什么样子：

​

我试图找出我做错了什么，但经过相当一段时间的测试，我回到了第一步。唯一让我担心的是，也许Matlab函数的工作方式与Mathcad中使用的函数有点不同(不能真正说明我完全理解它)。任何善良的灵魂都能告诉我它是否与fft函数有关？还是我只是个瞎子看不见他犯的一个愚蠢的错误？干杯，为我糟糕的英语道歉。

编辑

过了一段时间后，我的老师让我检查这种方法是否适用于相关的(某种)噪声，因为它同样适用于数学。相关后，其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中添加一些快速更改后的外观(一些微小的差异，但总的来说，它显示了我应该得到的效果)。它切断了扫描的乞求，但与我在这篇文章的开头所写的完全相同。

EDIT2

Nvm，它只是fft函数的一个维数问题。在将其转换为fft(dataw，[]，2)后，它看起来更好。

​

​

Answer 1

它不工作的主要原因是由于MathCad和MATLAB之间的快速傅立叶变换的比例因子。对于MathCad，有一个额外的缩放因子1/sqrt(N)，而MATLAB 不包括这个所谓的缩放因子。因此，如果您想要模拟使用MathCad看到的结果，则需要将快速傅立叶变换的结果乘以这个缩放因子。

此外，我对您的代码有一些建议：

1. 我们可以在没有任何循环的情况下完全实现这个矢量化。
2. 像fft和randn这样的函数可以对矩阵进行操作，您可以具体地将该函数应用于一个特定的维度。

请注意，我已经将您的随机噪声分布替换为Poisson随机噪声(来自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]);

因此，我们得到：

​

​