Преобразование Фурье с использованием 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))

Почему он разделяется на две вершины?

Answer 1

Вот несколько советов:

• использование np.fft.fft
• FFT начинается с 0 Гц
• нормализуют / отмасштабировать

Пример:

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

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

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]

x1 = 0.0
x2 = 5.0
N = 10000
T = x2 - x1

x = np.linspace(x1, x2, N)
y = light_intensity(x, 10.0, 0.0)
xf, yf = norm_fft(y, T, T / np.pi)

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 light_intensity(x, signal_gain, noise_gain):
    signal = norm.pdf(x, 0, 1)
    noise = np.random.randn(x.size)
    return signal_gain * signal + noise_gain * noise

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]

x1 = -10.0
x2 = 10.0
N = 10000
T = x2 - x1

x = np.linspace(x1, x2, N)
y = light_intensity(x, 10.0, 0.0)
xf, yf = norm_sym_fft(y, T, 4 / np.pi)

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

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

Это результат: