Python绘制fft.rfft频率的频率图

Question

这是我在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相同。但是x轴没有显示频率。我已经尝试了一些不同的方法，但都不能解决问题。 例如，在fft（这里我使用的是rfft）的情况下，应该可以解决这个问题。

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

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

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

Piero

- basetta81

Answer 1

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

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

- unutbu