计算任意信号的 DFT

Question

作为大学信号处理课程的一部分，我们被要求在 Matlab 中编写一个算法来使用 DFT 计算我们信号的单边频谱，而不使用内置的 fft() 函数matlab.这不是课程的评估部分，我只是想为自己获得这个 "right"。我目前使用的是 2018b 版的 Matlab，如果有人觉得这个有用的话。

我建立了一个 1 KHz 和 2KHz 正弦信号，相移 135 度（2*pi/3 弧度）。

然后使用 离散时间信号处理 (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 作为输入，在本例中为前面提到的信号，以及变换的分辨率，由 DFT_N 表示，已初始化为 100。此输出函数 X 应该是频域中的某个东西，但是绘制 X 会产生一个比单位圆略大的圆形图，并且在左手边有一个间隙。

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

非常感谢，J。

Answer 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-1 替换为 m=1:DFT_N 以表达担忧。

也用 2*pi*(n-1)*(m-1) 替换任何 2*pi*n*m 以获得正确的相位，因为 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));

Answer 2

如果您只对单边谱感兴趣，那么您可以计算 X(m) for m=1:DFT_N/2 因为 X 是围绕 m=DFT_N/2 共轭对称的，即 X(DFT_N/2+m) = X(DFT_N/2-m)'，由于 exp(-j*(pi*n+2*pi/DFT_N*m)) = exp(-j*(pi*n-2*pi/DFT_N*m))'。

作为旁注，对于给定的 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

计算任意信号的 DFT

Computing the DFT of an arbitrary signal

matlab

fft

frequency-analysis

dft