DFT 理解与python实现

我们知道傅里叶变换能够将时域的信号变换到频域分析，但是我们在计算机应用时，常常是离散的信号。如，用ADC采集回来的模拟信号，就变成了离散的信号。这时我们要对信号进行分析，就不能再使用傅里叶变换了，而要使用离散傅里叶变换，DFT：
$X[k]={\sum }_{n=0}^{N-1}x[n]e^{-j2\pi \ast kn/N}$
其中:

$X[k]$：离散频率下标为k时的频率大小
$x[n]$: 离散时域信号序列
$N$: 信号序列的长度，也就是采样的个数

下面我们用python来实现DFT：
对于信号$y=sin(2\ast \pi \ast f_0\ast t)$ 来进行DFT。
$f_0$ 设为10Hz
采样频率$f_s$设为50hz
采样个数N 设为100

所以$y[n]=x[n]=sin(2\ast \pi \ast f_0\ast n/f_s)$ $n=0\dots N-1$

$X[k]=sin(2\ast \pi \ast f_0\ast n/f_s)e^{-j2\pi \ast kn/N}$ $k=1\dots N-1$
code:

from numpy import arange, sin, pi, cos
import matplotlib.pyplot as plt
import cmath
from scipy.fftpack import fft,ifft

Kl = arange(0, 100, 1)
Xk = []
f0 = 10
fs = 50
N = len(Kl)
y = sin(2*pi*f0*Kl/fs)    # x[n]
for k in Kl:  # 求 X[k] k = 1....N
    temp = 0
    for l in range(0, N):  # 求X[k]
        cm = 1j
        cm *= 2*pi*l*k/N
        temp += sin(2*pi*f0*l/fs)*cmath.exp(cm)
    temp = temp
    Xk.append(abs(temp))


yy = fft(y)
yf = abs(yy)      # 用python自带的FFT验证
plt.plot(Kl, Xk, '*', label="DFT")
plt.hold(True)
plt.plot(Kl, yf, '>', label="python FFT")
plt.hold(True)
plt.plot(Kl, y, '', label="sin")
plt.legend(loc="lower right")
plt.show()

对比一下python自带FFT，发现自己计算出来的结果是正确的。

其中我们会发现在20这个点特别突出，是为什么呢？
我们会发现信号本身是10Hz的，采样频率是50Hz，所以信号是采样信号的1/5，而采样数为100，所以这里的20对应的频率就是10Hz。
另外在80这里也有一个突出，具体原因作者还没找到，暂时猜测是因为频谱泄露的原因，有人知道这是什么原因麻烦告诉一下作者好吗？ 后面待我弄清楚后，再更改这一部分。

(80突出是因为FFT对称性质)