Scipy FFT 和 Numpy FFT 在脉冲序列频谱上存在分歧？

Posted 2023-02-15

【中文标题】Scipy FFT 和 Numpy FFT 在脉冲序列频谱上存在分歧？

【英文标题】：Scipy FFT and Numpy FFT disagree on pulse train spectrum?

【发布时间】：2022-01-18 04:39:58

【问题描述】：

我正在对一系列脉冲进行 FFT。该系列是每 7 天一个幅度为 1 的脉冲，共 367 天。

当我运行以下代码时：

import numpy as np
import pandas as pd
from scipy.fft import fft, fftfreq, fftshift, ifft
from scipy.signal import blackman
from matplotlib import pyplot as plt
import random

## Signal 
num_samples = 367
# time in days
t = np.arange(int(num_samples))
# Amplitude and position of pulse. Amplitude here is 0 or 1 but can generate random values
# Position here is every 7th day
signal = [random.randint(1,1) if (i%7 == 0) else 0 for i, x in enumerate(t)]#np.sin(2*np.pi*5*t/N)#[random.randint(1,1) if (i%7 == 0) else 0 for i, x in enumerate(t)]#

# FFT and IFFT using Numpy

sr = 367
X = np.fft.fft(signal)
n = np.arange(num_samples)
T = num_samples/sr
freq = n/T 

plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.title('FFT using Numpy')
plt.stem(freq, np.abs(X), 'b', markerfmt=" ", basefmt="-b")
plt.xlabel('Freq (Hz)')
plt.ylabel('FFT Amplitude |X(freq)|')

plt.subplot(122)
plt.title('IFFT using Numpy')
plt.plot(t, np.fft.ifft(X), 'r')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()

# FFT and IFFT using Scipy

sp = fft(signal)
freq = fftfreq(t.shape[-1])

plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.title('FFT using Scipy')
plt.stem(freq, np.abs(sp), 'b', markerfmt=" ", basefmt="-b")
plt.xlabel('Freq (Hz)')
plt.ylabel('FFT Amplitude |sp(freq)|')

plt.subplot(122)
plt.title('IFFT using Scipy')
plt.plot(t, ifft(sp), 'r')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()

我得到以下信息：

显然存在移位和缩放问题，但更重要的是，我预计脉冲序列的 fft 是频谱中的一系列均匀峰值。我不明白导致的峰值，这意味着我可能误解了函数是如何解释信号的。任何指导将不胜感激。

【参考方案1】：

您应该使用 mp.fft.fftshift 来绘制 numpy fft。此外，还有一种更简单的方法可以使频率轴成为 numpy 情况。

你的代码用“>>“表示修改：

import numpy as np
import pandas as pd
from scipy.fft import fft, fftfreq, fftshift, ifft
from scipy.signal import blackman
from matplotlib import pyplot as plt
import random

## Signal 
num_samples = 367
# time in days
t = np.arange(int(num_samples))
# Amplitude and position of pulse. Amplitude here is 0 or 1 but can generate random values
# Position here is every 7th day
signal = [random.randint(1,1) if (i%7 == 0) else 0 for i, x in enumerate(t)]#np.sin(2*np.pi*5*t/N)#[random.randint(1,1) if (i%7 == 0) else 0 for i, x in enumerate(t)]#

# FFT and IFFT using Numpy

sr = 367

X = np.fft.fft(signal)  
freq = np.fft.fftfreq(len(t), d=1) # <<< Let numpy build the frequency axis for you >>>

plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.title('FFT using Numpy')

# <<< Shift the zero-frequency component to the center of the spectrum on frequ and X >>>
# <<< Note: use_line_collection=True in plt.stem() removes a warning (not important) >>>
plt.stem(np.fft.fftshift(freq), np.fft.fftshift(np.abs(X)), 'b', markerfmt=" ", basefmt="-b", use_line_collection=True)
plt.xlabel('Freq (Hz)')
plt.ylabel('FFT Amplitude |X(freq)|')

plt.subplot(122)
plt.title('IFFT using Numpy')
plt.plot(t, np.fft.ifft((X)), 'r') 
plt.plot(t, np.fft.ifft((X)), 'r') 

plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()

# FFT and IFFT using Scipy

sp = fft(signal)  

freq = fftfreq(len(t))
print(freq.shape)

plt.figure(figsize = (12, 6))
plt.subplot(121)
plt.title('FFT using Scipy')
plt.stem(freq, np.abs(sp), 'b', markerfmt=" ", basefmt="-b",use_line_collection=True)
plt.xlabel('Freq (Hz)')
plt.ylabel('FFT Amplitude |sp(freq)|')

plt.subplot(122)
plt.title('IFFT using Scipy')
plt.plot(t, ifft(sp), 'r')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()

输出：