fft.rfft 的 Python 绘图频率答案

【问题标题】：Python plot frequency of fft.rfftfft.rfft 的 Python 绘图频率
【发布时间】：2013-06-11 23:37:09
【问题描述】：

这是我在 stackoverflow 上的第一个问题，我希望我不会犯大错。我正在分析一组采样率为 1 Hz 的时间序列。为了研究它们的光谱，我需要绘制它们的傅里叶变换。

这是我的一段代码：

from obspy.core import read
import numpy as np 
import matplotlib.pyplot as plt

st = read('../SC_noise/*HEC_109C*_s', format='SAC')
stp = st.copy()
stp.detrend('linear')
stp.taper('cosine')

for tr in stp:
  dataonly = tr.data
  spec = np.fft.rfft(dataonly)
  plt.plot(abs(spec))
  plt.show()

这很好用：情节与我使用 SAC 得到的情节相同。但是 xaxis 不显示频率。我四处游荡了一下，发现了不同的想法：它们都不起作用。例如，在 fft 的情况下（这里我使用的是 rfft），这应该可以完成这项工作

samp_rate=1
freq = np.fft.fftfreq(len(spec), d=1./samp_rate)

但如果我使用它，它会给我带来负频率。

有人有想法吗？非常感谢您的所有帮助！

皮耶罗

【问题讨论】：

标签： python plot fft analysis frequency

【解决方案1】：

如果您的 NumPy 版本足够新（1.8 或更高），请使用 numpy.fft.rfftfreq。否则，这里是the definition：

def rfftfreq(n, d=1.0):
    """
Return the Discrete Fourier Transform sample frequencies
(for usage with rfft, irfft).

The returned float array `f` contains the frequency bin centers in cycles
per unit of the sample spacing (with zero at the start). For instance, if
the sample spacing is in seconds, then the frequency unit is cycles/second.

Given a window length `n` and a sample spacing `d`::

f = [0, 1, ..., n/2-1, n/2] / (d*n) if n is even
f = [0, 1, ..., (n-1)/2-1, (n-1)/2] / (d*n) if n is odd

Unlike `fftfreq` (but like `scipy.fftpack.rfftfreq`)
the Nyquist frequency component is considered to be positive.

Parameters
----------
n : int
Window length.
d : scalar, optional
Sample spacing (inverse of the sampling rate). Defaults to 1.

Returns
-------
f : ndarray
Array of length ``n//2 + 1`` containing the sample frequencies.

Examples
--------
>>> signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5, -3, 4], dtype=float)
>>> fourier = np.fft.rfft(signal)
>>> n = signal.size
>>> sample_rate = 100
>>> freq = np.fft.fftfreq(n, d=1./sample_rate)
>>> freq
array([ 0., 10., 20., 30., 40., -50., -40., -30., -20., -10.])
>>> freq = np.fft.rfftfreq(n, d=1./sample_rate)
>>> freq
array([ 0., 10., 20., 30., 40., 50.])

"""
    if not (isinstance(n,int) or isinstance(n, integer)):
        raise ValueError("n should be an integer")
    val = 1.0/(n*d)
    N = n//2 + 1
    results = arange(0, N, dtype=int)
    return results * val

【讨论】：