计算任意信号的 DFT答案

【问题标题】：Computing the DFT of an arbitrary signal计算任意信号的 DFT
【发布时间】：2018-10-22 19:58:30
【问题描述】：

作为大学信号处理课程的一部分，我们被要求在 Matlab 中编写一个算法，以使用 DFT 计算我们信号的单面频谱，而不使用 matlab 内置的 fft() 函数。这不是课程的评估部分，我只是有兴趣让自己“正确”。我目前正在使用 2018b 版本的 Matlab，如果有人觉得这很有用。

我构建了一个 1 KHz 和 2KHz 正弦波信号，相移 135 度 (2*pi/3 rad)。

然后使用离散时间信号处理（Allan V. Oppenheim）的9.1中的方程和欧拉公式来简化指数，我产生了这个代码：

%%DFT(currently buggy)
n=0;m=0;
for m=1:DFT_N-1 %DFT_Fmin;DFT_Fmax; %scrolls through DFT m values (K in text.)
    for n=1:DFT_N-1;%;(DFT_N-1);%<<redundant code? from Oppenheim eqn. 9.1 % eulers identity, K=m and n=n
        X(m)=x(n)*(cos((2*pi*n*m)/DFT_N)-j*sin((2*pi*n*m)/DFT_N));
        n=n+1;
    end
    %m=m+1; %redundant code?
end

这将 x 作为输入，在这种情况下是前面提到的信号，以及变换的分辨率，由已初始化为 100 的 DFT_N 表示。此函数的输出 X 应该是在频域中的某些东西，但绘制 X 会产生一个比单位圆稍大的圆形图，并且在左侧边缘有一个间隙。

我正在努力了解我应该如何将其转换为内置 DFT 算法给出的 stem() 图。

非常感谢，J.

【问题讨论】：

X 通常会很复杂，因此是圆形图。如果您对频率（幅度）响应感兴趣，通常需要绘制绝对值。
据我了解，x 绝对值的茎图将给出双面 DFT，然后添加绘图约束以获得单面？编辑
请注意 (a) Matlab 已经知道如何计算复指数，您不必通过扩展欧拉恒等式使代码复杂化。 (b) Matlab 函数在向量上运行得很好，所以你不需要嵌套循环。
（不允许我编辑）这会在 x 轴上的所有 n 的 y 轴上生成一个全局 X=1.107 的图，本质上是一条平线。这让我觉得我的绘图技术有其他问题
谢谢 Ben，我通常使用 C，所以对我来说嵌套循环更直观，尽管我知道它们不是在 matlab 中预编译的，因此速度较慢。我们从来没有上过matlab本身的课程，所以从来没有去正确理解语法，所以我只是一点一点地学习......至于欧拉展开式，我考虑过使用实部和虚部对于其他一些事情，但我再次看到它是多余的，因为 matlab 已经有一个系统来分离实部和虚部。

标签： matlab fft dft frequency-analysis

【解决方案1】：

这是你的错误：

将X(m)=x(n)*(cos.. 替换为X(m)=X(m)+x(n)*(cos..

对于给定的 m，它不会对变量 n 进行积分，而是仅覆盖 X(m)，仅覆盖 n = DFT_N-1 的最后一次计算。

请注意，对n=1:DFT_N-1 进行积分会省略一个谐波，即第一个谐波 exp(-j*2*pi)。代替 n=1:DFT_N-1 和 n=1:DFT_N 包括在内。我也会用 m=1:DFT_N 替换 m=1:DFT_N-1 来绘制问题。

还可以将任何2*pi*n*m 替换为2*pi*(n-1)*(m-1) 以获得正确的相位，因为X 的第一个条目应该对应于零频率，产生 sum_n x(n) * (cos(0) + j sin(0) ) = sum_n x(n)。如果您的信号 x 是实值，则零频率分量 X(1) 应该是实值，angle(X(1))=0。

最后一句话，不要忘记将零频率分量移到频谱中心以获得更好的可见性，X = circshift(X,floor(size(X)/2));

【讨论】：

【解决方案2】：

如果您只对单面光谱感兴趣，那么您可以只计算 X(m) for m=1:DFT_N/2，因为 X 在 m=DFT_N/2 周围是共轭对称的，即 X(DFT_N/2+m) = X(DFT_N/2-m)'，由于 @987654323 @。

附带说明，对于给定的 m，此程序计算数组 x 和另一个复指数数组之间的内积，即 exp(-j*2*pi/DFT_N*m*n)，对于 n = 0,1,...,N-1。 MATLAB 语法对于此类计算非常方便，您可以通过以下命令避免这种内循环

exp(-j*2*pi/DFT_N*m*(0:DFT_N-1)) * x

其中 x 是列向量。同样，您也可以通过为每个 m 逐行扩展复指数向量来避免第一个循环，即构建矩阵 exp(-j*2*pi/DFT_N*(0:DFT_N-1)'*(0:DFT_N-1))。那么你的 DFT 就是

X = exp(-j*2*pi/DFT_N*(0:DFT_N-1)'*(0:DFT_N-1)) * x

对于单面光谱，改为使用

X = exp(-j*2*pi/DFT_N*(0:floor((DFT_N-1)/2))'*(0:DFT_N-1)) * x

【讨论】：