Python plot frequency of fft.rfft

Question

this is my first question here on stackoverflow and I hope I will not make huge mistakes. I am analyzing a set of time series with sampling rate of 1 Hz. I need to plot their fourier transform in order to study their spectra.

Here it is my piece of code:

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()

This works just fine: the plot is the same I get using SAC. But the xaxis does not show frequencies. I've wandered around a little bit and found different ideas: none of them is working. For example in the case of a fft (here I am using a rfft) this should do the job

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

But if I use it it would give me negative frequencies.

Does anybody have an idea? Thank you very much in advance for all the help!

Piero

Answer 1

If your NumPy version is new enough (1.8 or better), use numpy.fft.rfftfreq. Otherwise, here is 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