使用Numpy进行傅里叶变换

Question

我正在尝试计算以下高斯函数的傅里叶变换：

# sample spacing
dx = 1.0 / 1000.0

# Points
x1 = -5
x2 = 5

x = np.arange(x1, x2, dx)

def light_intensity():
    return 10*sp.stats.norm.pdf(x, 0, 1)+0.1*np.random.randn(x.size)

fig, ax = plt.subplots()
ax.plot(x,light_intensity())

我在空间频域中创建了一个新的数组（高斯傅里叶变换是高斯分布，因此这些值应该是相似的）。我绘制出来，得到了以下结果：

fig, ax = plt.subplots()

xf = np.arange(x1,x2,dx)
yf= np.fft.fftshift(light_intensity())
ax.plot(xf,np.abs(yf))

为什么会分裂成两个峰值？

- Luke Polson

---

好的，这花了我一些时间，但基本上，你的x标签应该从0开始。 - Mateen Ulhaq

fftshift函数并不执行傅里叶变换。若要计算FFT，请使用numpy.fft.fft函数。fftshift函数可以对FFT结果进行移位操作。 - Warren Weckesser

当我将fft.fft替换到上面的代码中时，得到的图形更加奇怪。我在每一侧得到了两个delta函数峰值。 - Luke Polson

Answer 1

建议：

• 使用np.fft.fft
• FFT从0 Hz开始
• 对数据进行归一化/重新缩放

完整示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

def norm_fft(y, T, max_freq=None):
    N = y.shape[0]
    Nf = N // 2 if max_freq is None else int(max_freq * T)
    xf = np.linspace(0.0, 0.5 * N / T, N // 2)
    yf = 2.0 / N * np.fft.fft(y)
    return xf[:Nf], yf[:Nf]

def generate_signal(x, signal_gain=10.0, noise_gain=0.0):
    signal = norm.pdf(x, 0, 1)
    noise = np.random.randn(x.size)
    return signal_gain * signal + noise_gain * noise

# Signal parameters
x1 = 0.0
x2 = 5.0
N = 10000
T = x2 - x1

# Generate signal data
x = np.linspace(x1, x2, N)
y = generate_signal(x)

# Apply FFT
xf, yf = norm_fft(y, T, T / np.pi)

# Plot
fig, ax = plt.subplots(2)
ax[0].plot(x, y)
ax[1].plot(xf, np.abs(yf))
plt.show()

或者，带有噪声的情况:

---

对称图形

或者，如果您想要在频域中欣赏对称性:

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

def norm_sym_fft(y, T, max_freq=None):
    N = y.shape[0]
    b = N if max_freq is None else int(max_freq * T + N // 2)
    a = N - b
    xf = np.linspace(-0.5 * N / T, 0.5 * N / T, N)
    yf = 2.0 / N * np.fft.fftshift(np.fft.fft(y))
    return xf[a:b], yf[a:b]

def generate_signal(x, signal_gain=10.0, noise_gain=0.0):
    signal = norm.pdf(x, 0, 1)
    noise = np.random.randn(x.size)
    return signal_gain * signal + noise_gain * noise

# Signal parameters
x1 = -10.0
x2 = 10.0
N = 10000
T = x2 - x1

# Generate signal data
x = np.linspace(x1, x2, N)
y = generate_signal(x)

# Apply FFT
xf, yf = norm_sym_fft(y, T, 4 / np.pi)

# Plot
fig, ax = plt.subplots(2)
ax[0].plot(x, y)
ax[1].plot(xf, np.abs(yf))
plt.show()

或者，带有噪音的：

- Mateen Ulhaq

---

太棒了！感谢你的帮助！ - Luke Polson

快速问题。高斯函数的傅里叶变换应该是相同的高斯函数本身，对吗？例如：http://www.wolframalpha.com/input/?i=fourier+transform+of+4e%5E(-x%5E2%2F2)，其中傅里叶变换函数与原始函数（在频域中）相同。你展示的图中，fft函数的高度和标准差不同，这是有原因的吗？ - Luke Polson

我相信这是因为FFT=DFT，而不同于FT。 - Mateen Ulhaq

Answer 2

首先使用np.fft.fft来计算傅里叶变换，然后使用np.fft.fftshift将零频率分量移到频谱的中心。

用以下代码替换您代码的第二部分:

xf = np.arange(x1,x2,dx)
yf = np.fft.fft(light_intensity())
yfft = np.fft.fftshift(np.abs(yf))
fig,ax = plt.subplots(1,2,figsize=(10,5))
ax[0].plot(xf,light_intensity())
ax[1].plot(xf,yfft)
ax[1].set_xlim(-0.05,0.05)
plt.show()

这是结果：

- Mahdi