制作一个函数来进行傅里叶逆变换
Making a function to take the inverse fourier transform
我正在尝试通过制作自己的函数来进行傅里叶逆变换。这是对我的时间序列进行傅里叶变换的函数,它看起来工作正常。
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
我的功能有什么问题?我没有收到任何错误,但返回的信号完全错误。
在你的倒数中,指数应该是正的(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()
给出(省略第一个正弦图):
和
运行 你编码你应该得到以下警告:
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
函数)。
最后,如GrimTrigger :
your inverse the exponent should be positive (2j instead of -2j) and drop the /N
我正在尝试通过制作自己的函数来进行傅里叶逆变换。这是对我的时间序列进行傅里叶变换的函数,它看起来工作正常。
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
我的功能有什么问题?我没有收到任何错误,但返回的信号完全错误。
在你的倒数中,指数应该是正的(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()
给出(省略第一个正弦图):
和
运行 你编码你应该得到以下警告:
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
函数)。
最后,如GrimTrigger
your inverse the exponent should be positive (2j instead of -2j) and drop the /N