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を使用して取得同じです。しかし、xaxisは周波数を表示しません。私は少し歩き回り、さまざまなアイデアを見つけました：誰も働いていません。 FFTの場合、例えば （ここで私はRFFTを使用しています）これは

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

仕事をする必要があります。しかし、私はそれを使用する場合、それは私に負の周波数を与えるだろう。

誰かがアイデアを持っていますか？ 大変ありがとうございました！

ピエロ

Answer 1

numpy.fft.rfftfreqを使用して、あなたのnumpyのバージョンは、（1.8以上）に十分な新規である場合。そうでない場合は、ここに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