dft算法理解及实现

前面对傅里叶变化的含义有了一定的理解，现在，就自己来实现一下DFT算法，并记录下遇到的问题。

一、DFT算法简述

DFT算法是一种最基础的离散傅里叶变换，其计算公式为：
{ x ( k ) } = ∑ n = 0 N − 1 x ( n ) ⋅ e − j 2 π N n k , k = 0 , 1 , 2 , . . . , N − 1 \{x(k)\}=\sum_{n=0}^{N-1}x(n)\cdot e^{-j\frac{2\pi}{N}nk}, k=0, 1, 2,...,N-1 {x(k)}=n=0∑N−1​x(n)⋅e−jN2π​nk,k=0,1,2,...,N−1
其中，k即为旋转因子$W_n$的旋转系数。也就是说，它是直接使用$w=2\pi k$的正弦函数和余弦函数进行逐个尝试的。

二、代码实现

首先是定义dft函数

% 傅里叶变换，用于将时域数据转为频域数据

% 输入为：时域数据、频域采样点数量
function freq_domain = dft(time_domain, freq_N)

time_N = length(time_domain);
freq_domain = zeros(1, freq_N);
x = 0:time_N-1;
x = x / time_N;

for k=0:freq_N-1
    Wn = exp(-1i * 2 * pi * x * k);
    freq_domain(k + 1) = Wn * time_domain';
end

end

然后直接调用，并于fft()函数的结果进行比较：

clc, clear

time_N = 200;
x = 0:time_N-1;
x = x / time_N;

y = 3 * sin(2 * pi * 30 * x) + 4 * sin(2 * pi * 60 * x);

w1 = dft(y, time_N);
w2 = fft(y);

w11 = abs(w1);
w11(w11<1)=0;
f1 = find(w11)
w_dft = w11(f1) / time_N * 2

w22 = abs(w2);
w22(w22<1)=0;
f2 = find(w22)
w_fft = w22(f1) / time_N * 2

幅值和相位都与理论的一致。

三、存在问题

但是，当查看w1和w2的值时，会发现它们其实并不一致。

当把调用的代码中，x的时间归一化进行删除时，结果且基本一致，但此时，已经对相位进行了修改，无法在计算出所对应的幅值了，所以意义也不大。

clc, clear

time_N = 200;
x = 0:time_N-1;
% x = x / time_N;

y = 3 * sin(2 * pi * 30 * x) + 4 * sin(2 * pi * 60 * x);

同时，查看matlab的fft代码解释，依旧没有找出问题所在。

%FFT Discrete Fourier transform.
%   FFT(X) is the discrete Fourier transform (DFT) of vector X.  For
%   matrices, the FFT operation is applied to each column. For N-D
%   arrays, the FFT operation operates on the first non-singleton
%   dimension.
%
%   FFT(X,N) is the N-point FFT, padded with zeros if X has less
%   than N points and truncated if it has more.
%
%   FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the
%   dimension DIM.
%   
%   For length N input vector x, the DFT is a length N vector X,
%   with elements
%                    N
%      X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
%                   n=1
%   The inverse DFT (computed by IFFT) is given by
%                    N
%      x(n) = (1/N) sum  X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
%                   k=1
%
%   See also FFT2, FFTN, FFTSHIFT, FFTW, IFFT, IFFT2, IFFTN.

%   Copyright 1984-2005 The MathWorks, Inc.

%   Built-in function.

四、总结

目前依旧无法定位到问题所在，或许是我哪里理解有问题，有空接着搞了。