【问题标题】:Plotting Fourier Transform of Gaussian function with python, but the result was wrong用python绘制高斯函数的傅里叶变换,但结果错误
【发布时间】:2020-01-30 10:24:57
【问题描述】:

我想了很久,但我不知道问题是什么。希望你能帮助我,谢谢。

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()

结果:

【问题讨论】:

  • 你期望它是什么?
  • 我想得到高斯函数傅里叶变换后的曲线
  • 与您得到的答案相比,这应该如何?
  • 高斯函数在 FT 之后应该仍然是高斯函数
  • 你不应该取结果的实部,而是它的绝对值

标签: python fft gaussian


【解决方案1】:

如前所述,您需要绝对值,而不是实数部分。 一个最小的例子,显示 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');

【讨论】:

  • 非常感谢您的帮助
【解决方案2】:

你必须从时间尺度改变到频率尺度

【讨论】:

【解决方案3】:

当您进行 FFT 时,您将获得 simetric 变换,即正曲线到负曲线的镜像。通常,您只会看到积极的一面。 此外,您应该注意采样率,因为 FFT 旨在将 时域 输入转换为 频域,输入信息的时间或采样率很重要。因此,在np.fft.fftfreq(n, d=timestep) 中为您的采样率添加时间步长。

如果您只是想制作一个 fft 的正常 dist 信号,这里有另一个问题,以及关于您为什么会出现这种行为的一些很好的解释:

Fourier transform of a Gaussian is not a Gaussian, but thats wrong! - Python

【讨论】:

    【解决方案4】:

    你的代码有两个错误:

    • 绘图时不要取实数,取绝对值。
    • 来自文档:

    如果A = fft(a, n),则A[0] 包含零频率项(均值 信号),对于实际输入,它总是纯真实的。然后 A[1:n/2] 包含正频率项,A[n/2+1:] 包含 负频率项,按负数递减的顺序 频率。

    您可以使用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 轴。

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2014-06-10
      • 1970-01-01
      • 1970-01-01
      • 2014-11-02
      • 1970-01-01
      相关资源
      最近更新 更多