【问题标题】:Scipy FFT and Numpy FFT disagree on pulse train spectrum?Scipy FFT 和 Numpy FFT 在脉冲序列频谱上存在分歧?
【发布时间】: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 是频谱中的一系列均匀峰值。我不明白导致的峰值,这意味着我可能误解了函数是如何解释信号的。任何指导将不胜感激。

【问题讨论】:

    标签: python-3.x numpy scipy 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()
    

    输出:

    【讨论】:

      猜你喜欢
      • 2019-01-24
      • 1970-01-01
      • 2020-10-23
      • 2013-03-12
      • 1970-01-01
      • 2012-04-30
      • 2021-03-09
      • 1970-01-01
      • 2015-07-16
      相关资源
      最近更新 更多