Advertisement

离散傅里叶变换DFT

阅读量:

离散傅里叶变换DFT


DFT和IDFT

时域中长度为N的序列x[n]的离散傅里叶变换(DFT)和逆变换(IDFT)

为了方便记忆,常用一个中间变量

这样DFT和IDFT变成下面容易记忆的形式:

进一步,正反变换都可以看成是两个矩阵的乘法

观察到对 W 矩阵进行求逆操作转换为 W^{-1} 的过程只需将每个元素取倒数,并除以 即可完成。无需直接计算其逆矩阵,在理论上这种转换过程非常高效。


numpy实现

代码参考自这里

波形显示函数

复制代码
    def show(ori_func, ft, sampling_period = 5): 
    n = len(ori_func) 
    interval = sampling_period / float(n) 
    print interval
    # 绘制原始函数
    plt.subplot(2, 1, 1) 
    plt.plot(np.arange(0, sampling_period, interval), ori_func, 'black') 
    plt.xlabel('Time'), plt.ylabel('Amplitude') 
    # 绘制变换后的函数
    plt.subplot(2,1,2) 
    frequency = np.arange(n / 2) / (n * interval) 
    nfft = abs(ft[range(int(n / 2))] / n ) 
    plt.plot(frequency, nfft, 'red') 
    plt.xlabel('Freq (Hz)'), plt.ylabel('Amp. Spectrum') 
    plt.show() 

多个谐波叠加波形

复制代码
    time = np.arange(0, 5, 0.005) 
    x1 = np.sin(2 * np.pi * 1 * time) # 频率为1
    x2 = np.sin(2 * np.pi * 20 * time) # 频率为20
    x3 = np.sin(2 * np.pi * 60 * time) # 频率为60
    x = x1 + x2 + x3 # 叠加
    y = np.fft.fft(x) 
    show(x, y) 
这里写图片描述

验证numpy.fft.fft的计算原理

复制代码
    # 傅里叶变换DFT
    x = np.random.random(500) 
    N = len(x) 
    n = np.arange(N) 
    k = n.reshape((N, 1)) 
    W = np.exp(-2j * np.pi * k * n / N) # 500*500 矩阵
    y = np.dot(W, x) 
    print np.allclose(y, np.fft.fft(x)) # 应该是 True

验证numpy.fft.ifft的计算原理

复制代码
    # 傅里叶逆变换IDFT
    W_inv = np.exp(2j * np.pi * k * n / N) / N
    x = np.dot(y, W_inv) 
    print np.allclose(x, np.fft.ifft(y)) # 应该是 True

参考资料

博客

博客

API

书籍

全部评论 (0)

还没有任何评论哟~