制作一个函数来进行傅里叶逆变换

Posted

技术标签:

【中文标题】制作一个函数来进行傅里叶逆变换【英文标题】:Making a function to take the inverse fourier transform 【发布时间】:2020-12-19 10:58:29 【问题描述】:

我正在尝试通过创建自己的函数来进行傅立叶逆变换。这是对我的时间序列进行傅里叶变换的函数,它看起来工作正常。

def DFT(x, frequencies):
    N1 = x.size
    k = frequencies.size
    t = np.arange(N1)
    fourier = np.zeros(k)
    for i in range(0,k):
        fourier[i] = np.dot(x,(1/k)*np.exp(-2j*np.pi*frequencies[i]*t))
    return fourier

这是我的原始信号(只是一个正弦波):

N2 = 1*10**6
time = np.arange(0,N2,1000)
lam = .1*10**6
f = 1/lam
freq = np.arange(0,.05,.0001)


signal = np.sin(2*np.pi*f*time)

功率谱是使用我的 DFT(傅立叶函数)绘制的:

plt.plot(freq, np.abs(DFT(signal,freq))**2)
plt.xlabel('Frequency')
plt.title('Spectrum of Sine Wave')
plt.grid()
plt.show()

但是当我尝试将我的函数应用于傅里叶逆变换时,我没有得到我原来的正弦波:

def IFT(fft, frequencies):
    N = fft.size
    k = frequencies.size
    n = np.arange(N)
    inverse_fourier = np.zeros(k)
    for i in range(0,k):
        inverse_fourier[i] = np.dot(fft,np.exp((-2j*np.pi*frequencies[i]*n)/N)) #[None,:]
    return inverse_fourier

我的功能有什么问题?我没有收到任何错误,但返回的信号完全错误。

【问题讨论】:

【参考方案1】:

运行您的代码,您应该会收到以下警告:

ComplexWarning: Casting complex values to real discards the imaginary part
  fourier[i] = np.dot(x,(1/k)*np.exp(-2j*np.pi*frequencies[i]*t))

由于生成的傅立叶变换应该是复数值,因此该警告应引起关注。要消除此警告,您可以像这样初始化fourier

fourier = np.zeros(k, dtype=complex)

Discrete Fourier Transform 的公式还包括覆盖整个[0,1) 范围的频率的总和。要获得 1000 点 DFT(正如您在代码中所做的那样),您必须使用

freq = np.arange(0,1,.001)

这将产生一个包含 2 个尖峰的频谱:一个处于预期频率,另一个对称的高于奈奎斯特频率。在绘制实值信号的频谱时,通常会丢弃高于奈奎斯特频率的结果(但在 IFT 函数中使用完整频谱)。

最后,作为GrimTriggerpointed out:

您的指数倒数应该是正数(2j 而不是 -2j)并删除 /N

【讨论】:

【参考方案2】:

在你的倒数中,指数应该是正数(2j 而不是-2j)并删除/N,它给出(添加了用于演示的图):

import numpy as np
import matplotlib.pyplot as plt

def DFT(x, frequencies):
    N1 = x.size
    k = frequencies.size
    t = np.arange(N1)
    fourier = np.zeros(k)
    for i in range(0,k):
        fourier[i] = np.dot(x, (1/k)*np.exp(-2j*np.pi*frequencies[i]*t))
    return fourier

def IFT(fft, frequencies):
    N = fft.size
    k = frequencies.size
    n = np.arange(N)
    inverse_fourier = np.zeros(k)
    for i in range(0,k):
        inverse_fourier[i] = np.dot(fft, np.exp((2j*np.pi*frequencies[i]*n))) #[None,:]
    return inverse_fourier

N2 = 1*10**6
time = np.arange(0,N2,2000)
lam = .1*10**6
f = 1/lam
freq = np.arange(0,.05,.0001)

signal = np.sin(2*np.pi*f*time)

plt.plot(time, signal)
plt.xlabel('Time')
plt.title('Sine Wave')
plt.grid()
plt.show()

dft = DFT(signal, freq)

plt.plot(freq, np.abs(dft)**2)
plt.xlabel('Frequency')
plt.title('Spectrum of Sine Wave')
plt.grid()
plt.show()

plt.plot(time, IFT(dft, freq))
plt.xlabel('Time')
plt.title('Sine Wave')
plt.grid()
plt.show()

给出(第一个罪图省略):

【讨论】:

以上是关于制作一个函数来进行傅里叶逆变换的主要内容,如果未能解决你的问题,请参考以下文章

2021-05-21 Matlab实现快速傅里叶逆变换

傅里叶逆变换FFT3W

音频算法入门-傅里叶变换

仅具有幅度的傅里叶逆变换 - 我也需要相位吗?

图像处理之图像傅里叶变换

傅里叶变换的结论汇总