问题描述
我正在尝试计算以下高斯的傅立叶变换：
# 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))
为什么会分裂成两个峰？
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()
或者，有噪音：
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()
这是结果：