7

これは、ここでのスタックオーバーフローに関する最初の質問であり、大きな間違いを犯さないことを願っています。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)

しかし、それを使用すると、負の周波数が得られます。

誰かアイデアがありますか？すべての助けを前もってありがとうございました！

ピエロ

4

1 に答える 1

6

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
于 2013-06-15T12:10:30.007 に答える