用python绘制高斯函数的傅立叶变换,但结果是错误的
Plotting Fourier Transform of Gaussian function with python, but the result was wrong
我想了很久,也没找到问题所在。希望你能帮助我,谢谢。
F(s) Gaussian function
F(s)=1/(√2π s) e^(-(w-μ)^2/(2s^2 ))
代码:
import numpy as np
from matplotlib import pyplot as plt
from math import pi
from scipy.fft import fft
def F_S(w, mu, sig):
return (np.exp(-np.power(w-mu, 2)/(2 * np.power(sig, 2))))/(np.power(2*pi, 0.5)*sig)
w=np.linspace(-5,5,100)
plt.plot(w, np.real(np.fft.fft(F_S(w, 0, 1))))
plt.show()
结果:
你必须从时间尺度改为频率尺度
当您进行 FFT 时,您将获得对称变换,即正曲线到负曲线的镜像。通常,你只会看到积极的一面。
此外,您应该注意采样率,因为 FFT 旨在将 时域 输入转换为 频域 、时间或采样率,输入信息很重要。因此,在 np.fft.fftfreq(n, d=timestep)
中为您的采样率添加时间步长。
如果你只是想对正常的 dist 信号进行 fft,这里有另一个问题和一些关于为什么你会出现这种行为的很好的解释:
Fourier transform of a Gaussian is not a Gaussian, but thats wrong! - Python
你的代码有两个错误:
- 不取实部,作图时取绝对值
- 来自文档:
If A = fft(a, n)
, then A[0]
contains the zero-frequency term (the mean
of the signal), which is always purely real for real inputs. Then
A[1:n/2]
contains the positive-frequency terms, and A[n/2+1:]
contains
the negative-frequency terms, in order of decreasingly negative
frequency.
您可以重新排列 np.fft.fftshift
的元素。
工作代码:
import numpy as np
from matplotlib import pyplot as plt
from math import pi
from scipy.fftpack import fft, fftshift
def F_S(w, mu, sig):
return (np.exp(-np.power(w-mu, 2)/(2 * np.power(sig, 2))))/(np.power(2*pi, 0.5)*sig)
w=np.linspace(-5,5,100)
plt.plot(w, fftshift(np.abs(np.fft.fft(F_S(w, 0, 1)))))
plt.show()
此外,您可能还想考虑缩放 x 轴。
如前所述,您需要的是绝对值,而不是实部。
一个最小的例子,显示 re/im , abs/phase 光谱。
import numpy as np
import matplotlib.pyplot as p
%matplotlib inline
n=1001 # add 1 to keep the interval a round number when using linspace
t = np.linspace(-5, 5, n ) # presumed to be time
dt=t[1]-t[0] # time resolution
print(f'sampling every {dt:.3f} sec , so at {1/dt:.1f} Sa/sec, max. freq will be {1/2/dt:.1f} Hz')
y = np.exp(-(t**2)/0.01) # signal in time
fr= np.fft.fftshift(np.fft.fftfreq(n, dt)) # shift helps with sorting the frequencies for better plotting
ft=np.fft.fftshift(np.fft.fft(y)) # fftshift only necessary for plotting in sequence
p.figure(figsize=(20,12))
p.subplot(231)
p.plot(t,y,'.-')
p.xlabel('time (secs)')
p.title('signal in time')
p.subplot(232)
p.plot(fr,np.abs(ft), '.-',lw=0.3)
p.xlabel('freq (Hz)')
p.title('spectrum, abs');
p.subplot(233)
p.plot(fr,np.real(ft), '.-',lw=0.3)
p.xlabel('freq (Hz)')
p.title('spectrum, real');
p.subplot(235)
p.plot(fr,np.angle(ft), '.-', lw=0.3)
p.xlabel('freq (Hz)')
p.title('spectrum, phase');
p.subplot(236)
p.plot(fr,np.imag(ft), '.-',lw=0.3)
p.xlabel('freq (Hz)')
p.title('spectrum, imag');
我想了很久,也没找到问题所在。希望你能帮助我,谢谢。
F(s) Gaussian function
F(s)=1/(√2π s) e^(-(w-μ)^2/(2s^2 ))
代码:
import numpy as np
from matplotlib import pyplot as plt
from math import pi
from scipy.fft import fft
def F_S(w, mu, sig):
return (np.exp(-np.power(w-mu, 2)/(2 * np.power(sig, 2))))/(np.power(2*pi, 0.5)*sig)
w=np.linspace(-5,5,100)
plt.plot(w, np.real(np.fft.fft(F_S(w, 0, 1))))
plt.show()
结果:
你必须从时间尺度改为频率尺度
当您进行 FFT 时,您将获得对称变换,即正曲线到负曲线的镜像。通常,你只会看到积极的一面。
此外,您应该注意采样率,因为 FFT 旨在将 时域 输入转换为 频域 、时间或采样率,输入信息很重要。因此,在 np.fft.fftfreq(n, d=timestep)
中为您的采样率添加时间步长。
如果你只是想对正常的 dist 信号进行 fft,这里有另一个问题和一些关于为什么你会出现这种行为的很好的解释:
Fourier transform of a Gaussian is not a Gaussian, but thats wrong! - Python
你的代码有两个错误:
- 不取实部,作图时取绝对值
- 来自文档:
If
A = fft(a, n)
, thenA[0]
contains the zero-frequency term (the mean of the signal), which is always purely real for real inputs. ThenA[1:n/2]
contains the positive-frequency terms, andA[n/2+1:]
contains the negative-frequency terms, in order of decreasingly negative frequency.
您可以重新排列 np.fft.fftshift
的元素。
工作代码:
import numpy as np
from matplotlib import pyplot as plt
from math import pi
from scipy.fftpack import fft, fftshift
def F_S(w, mu, sig):
return (np.exp(-np.power(w-mu, 2)/(2 * np.power(sig, 2))))/(np.power(2*pi, 0.5)*sig)
w=np.linspace(-5,5,100)
plt.plot(w, fftshift(np.abs(np.fft.fft(F_S(w, 0, 1)))))
plt.show()
此外,您可能还想考虑缩放 x 轴。
如前所述,您需要的是绝对值,而不是实部。 一个最小的例子,显示 re/im , abs/phase 光谱。
import numpy as np
import matplotlib.pyplot as p
%matplotlib inline
n=1001 # add 1 to keep the interval a round number when using linspace
t = np.linspace(-5, 5, n ) # presumed to be time
dt=t[1]-t[0] # time resolution
print(f'sampling every {dt:.3f} sec , so at {1/dt:.1f} Sa/sec, max. freq will be {1/2/dt:.1f} Hz')
y = np.exp(-(t**2)/0.01) # signal in time
fr= np.fft.fftshift(np.fft.fftfreq(n, dt)) # shift helps with sorting the frequencies for better plotting
ft=np.fft.fftshift(np.fft.fft(y)) # fftshift only necessary for plotting in sequence
p.figure(figsize=(20,12))
p.subplot(231)
p.plot(t,y,'.-')
p.xlabel('time (secs)')
p.title('signal in time')
p.subplot(232)
p.plot(fr,np.abs(ft), '.-',lw=0.3)
p.xlabel('freq (Hz)')
p.title('spectrum, abs');
p.subplot(233)
p.plot(fr,np.real(ft), '.-',lw=0.3)
p.xlabel('freq (Hz)')
p.title('spectrum, real');
p.subplot(235)
p.plot(fr,np.angle(ft), '.-', lw=0.3)
p.xlabel('freq (Hz)')
p.title('spectrum, phase');
p.subplot(236)
p.plot(fr,np.imag(ft), '.-',lw=0.3)
p.xlabel('freq (Hz)')
p.title('spectrum, imag');