在MATLAB中使用DFT产生的意外结果

问题描述：

在MATLAB中计算离散傅里叶变换时出现问题，显然得到正确的结果，但当绘制所获得频率的幅度时，您可以看到非常接近零的值，它应该完全为零。我用我自己的实现：我知道使用MATLAB的fft更好，但我需要使用我自己的实现，因为它是为大学。在MATLAB中使用DFT产生的意外结果

我用来产生方波代码：

x = [ones(1,8), -ones(1,8)]; 
for i=1:63 
    x = [x, ones(1,8), -ones(1,8)]; 
end

MATLAB版本：R2013a（8.1.0.604）64位

我已尝试发生在我身上的一切，但我没有很多使用MATLAB的经验，我还没有在论坛中发现与此问题有关的信息。我希望有一个人可以帮助我。

在此先感谢。

我认为这是一个数值问题。由于这些值在'1e-15'的范围内，而峰值的幅度大约为'656'，所以我不会打扰太多。 MATLAB FFT和你的DFT程序之间的平方误差的总和在'1e-20'的范围内，所以基本上为零。但也许别人有更详细的消除？ – hbaderts

我同意@hbaderts。这些值非常小，几乎和* eps *一样小，并且在这里可以被认为是0。你的FFT解决方案看起来很棒 – siliconwafer

感谢您的回答，我知道这个值是微不足道的，并不显示在情节，但是当我想绘制角度函数的相位获得无意义的结果。我使用y（（y -0.00001））= 0来隐藏这个值，但它不是一个优雅的解决方案。 –

答

这将是一个数值问题。值在1e-15的范围内，而信号的DFT的值在1e+02的范围内。进行进一步处理时，这很可能不会导致任何错误。您可以通过

y = fft(x); 
yh = Discrete_Fourier_Transform(x); 
sum(abs(yh - y).^2) 
ans = 
    3.1327e-20

这基本上是零计算你的DFT和MATLAB fft功能之间的总平方误差。因此我会得出结论：你的DFT函数工作得很好。

只是一个小小的评论：你可以很容易地矢量化DFT。

n = 0:1:N-1; 
k = 0:1:N-1; 
y = exp(-1j*2*pi/N * n'*k) * x(:);

随着n'*k您创建的n和k所有组合的矩阵。然后你采取这些矩阵元素中的每一个的exp(...)。使用x(:)确保x是一个列向量，因此您可以执行矩阵乘法运算(...)*x，它会自动汇总所有k的结果。其实，我只是注意到，这正是DFT众所周知的矩阵形式。