python - 実験データの周波数領域へのフーリエ変換

Question

私はpythonがまったく初めてです。波動データのテスト実験を行いました。時系列データを入手できます。周波数ドメインでそれを示すにはどうすればよいですか? 参照できる例はありますか？以下のようなプログラムを考えましたが、うまくいきません。助けてください。

#Program for Fourier Transformation
import numpy as np
import numpy.fft as fft
import matplotlib.pyplot as plt

def readdat( filename ):
    """
        Reads sectional area curve data from file filename
    """

    # read all lines of input files
    fp = open( filename, 'r')
    lines = fp.readlines() # to read the tabulated data
    fp.close()

    # interpret data
    time = []
    ampl = []
    for line in lines:
        if line[0:1] == '#':
            continue # ignore comments in the file
        try:
            time.append(float(line.split()[0])) #first column is time
            ampl.append(float(line.split()[1])) # second column is corresponding amplitude
        except:
            # if the data interpretation fails..
            continue
    return np.asarray(time), np.asarray(ampl)

if __name__ == '__main__':

    time, ampl = readdat( 'data.dat')
    print time
    print ampl

spectrum = fft.fft(ampl)
freq = fft.fftfreq(len(spectrum))
print freq

Answer 1

いくつかの結果をプロットするためのプログラムの最小限の修正は、次のようなものです

#Program for Fourier Transformation
import numpy as np
import numpy.fft as fft
import matplotlib.pyplot as plt

def readdat( filename ):
    """
        Reads sectional area curve data from file filename
    """

    # read all lines of input files
    fp = open( filename, 'r')
    lines = fp.readlines() # to read the tabulated data
    fp.close()

    # interpret data
    time = []
    ampl = []
    for line in lines:
        if line[0:1] == '#':
            continue # ignore comments in the file
        try:
            time.append(float(line.split()[0])) #first column is time
            ampl.append(float(line.split()[1])) # second column is corresponding amplitude
        except:
            # if the data interpretation fails..
            continue
    return np.asarray(time), np.asarray(ampl)

if __name__ == '__main__':

    time, ampl = readdat( 'data.dat')
    print time
    print ampl

spectrum = fft.fft(ampl)
timestep = time[1]-time[0] # assume samples at regular intervals
freq = fft.fftfreq(len(spectrum),d=timestep)
freq=fft.fftshift(freq)
spectrum = fft.fftshift(spectrum)
plt.figure(0,figsize=(5.0*1.21,5.0))
plt.plot(freq,spectrum)
print freq
plt.xlabel("frequencies")
plt.ylabel("spectrum")
plt.savefig("/tmp/figure.png")