使用python实现离散时间傅里叶变换

3,286 阅读3分钟

我们经常使用傅里叶变换来计算数字信号的频谱,进而分析数字信号,离散时间傅里叶变换的公式为:

X(e^{jw})=\sum_{n=-\infty}^{n=+\infty} x[n]e^{-jwn}

可是自己动手实现一遍才是最好的学习。
在数字分析里面,傅里叶变换默认等时间间隔采样,不需要时间序列,只需要信号数组即可分析。
分析过程如下:

  1. 对于含有 n 个样本值的数字信号序列,根据奈奎斯特采样定律,包含的周期数最大为n/2,周期数为 0 代表直流分量。所以,当周期数表示为离散的0,1,2,3...n/2,总的数目为n/2+1
  2. 傅里叶变换之后的结果为复数, 下标为k的复数a+b*j表示时域信号中周期为N/k个取样值的正弦波和余弦波的成分的多少, 其中a表示cos波形的成分,b表示sin波形的成分
  3. 首先产生一个长度为n,一倍周期的e^{-jwn}(即为cos(wn)-jsin(wn))波样本序列.
  4. 将数字信号序列中的每一个样本与1倍周期的样本波形序列相乘,得到n个乘积,将n个乘积相加,放入f[1]中。
  5. 再产生一个长度为n,两倍周期的e^{-jwn}(即为cos(wn)-jsin(wn))波样本序列,再将数字信号序列中的每一个样本与2倍周期的样本波形序列相乘,得到n个乘积,将n个乘积相加,放入f[2]中。依次重复。
  6. 对于0倍周期,即直流分量来说,可以认为产生的是0倍周期的样本波形,重复操作,放入f[0]即可。
  7. 这样就得到了数字信号序列的傅里叶变换

使用方法:
从以上过程得到数字序列的傅里叶变换之后,如果想要得到真正频谱,还需要做处理:

  1. 计算出的每一个频率下的幅值需要除以时间序列的长度,类似求平均的过程
  2. 每一个频率下的幅值是一个复数,需要对它求模,而且因为在负频率处也有值,所以需要对于实信号需要乘2
  3. 频率的序列为0到采样率的一半,长度为n/2+1 完整程序:
# 离散时间傅里叶变换的python实现
import numpy as np
import math
import pylab as pl
import scipy.signal as signal
import matplotlib.pyplot as plt

sampling_rate=1000
t1=np.arange(0, 10.0, 1.0/sampling_rate)
x1 =np.sin(15*np.pi*t1)

# 傅里叶变换
def fft1(xx):
#     t=np.arange(0, s)
    t=np.linspace(0, 1.0, len(xx))
    f = np.arange(len(xx)/2+1, dtype=complex)
    for index in range(len(f)):
        f[index]=complex(np.sum(np.cos(2*np.pi*index*t)*xx), -np.sum(np.sin(2*np.pi*index*t)*xx))
    return f

# len(x1)
xf=fft1(x1)/len(x1)
freqs = np.linspace(0, sampling_rate/2, len(x1)/2+1)
plt.figure(figsize=(16,4))
plt.plot(freqs,2*np.abs(xf),'r--')

plt.xlabel("Frequency(Hz)")
plt.ylabel("Amplitude($m$)")
plt.title("Amplitude-Frequency curve")

plt.show()

png


plt.figure(figsize=(16,4))
plt.plot(freqs,2*np.abs(xf),'r--')

plt.xlabel("Frequency(Hz)")
plt.ylabel("Amplitude($m$)")
plt.title("Amplitude-Frequency curve")
plt.xlim(0,20)
plt.show()

png

此处实现的是传统的傅里叶变换,这种方法实际已经不用了,现在使用快速傅里叶变换,其实两种是等价的,但是快速傅里叶变换时间复杂度要小很多。