hidden no source

継ぎ足し使う

【Python】f特抽出

正しくいうと抽出なんかじゃなくてピーク値を抽出するっていうプログラム

#coding:utf-8
import wave
import numpy as np
import scipy.fftpack
from scipy import signal
from pylab import *
import glob, os
#----------------------------------

fs = 16384#wf.getframerate()  # サンプリング周波数
start = 0  # サンプリングする開始位置
N =16384 #65536    # FFTのサンプル数



def get_filelist():
    #------------------------------+
    #ファイル名取得-ソート
    #------------------------------+
    global files,name
    files = glob.glob('*.csv')
    name  = []
    for x in files:
        no,ext = os.path.splitext(x)
        name.append(str(no))#[name]へ取得したファイル名リスト[no]から取得して格納していく
    
    name.sort()
    print files
    for x in xrange( len(name) ):
        name[x] = str(name[x]) + '.csv'

def get_peak(amp):
    #------------------------------+
    #ピーク値の抽出を行う.
    #in:FFT後の音圧レベル情報
    #out:ピーク位置の音圧レベル,周波数値を抜き出したリスト(総データ数は対応)
    #-------------------------------+
    frq_result = []
    lvl_result = []
    for i in xrange(2,len(amp) -1):
        if amp[i] > amp[i+1] and amp[i] > amp[i-1]:
            lvl_result.append( amp[i] )
            frq_result.append( freqList[i])

    print(len(lvl_result),len(frq_result) )
    return (lvl_result, frq_result)

get_filelist()


all_lvl = []
all_frq = []
for file in xrange(len(files)):
    filename = str(files[file])
    A,x = np.genfromtxt(filename,delimiter = ",",skiprows = 2,unpack = True ) 
    
    print file, len(files)
    if __name__ == "__main__" :
        global freqList
        X = np.fft.fft(x[start:start+N])  # FFT
        
        freqList = np.fft.fftfreq(N, d=1.0/fs)  # 周波数軸の値を計算
        amplitudeSpectrum = [20*log10 (np.sqrt(c.real ** 2 + c.imag ** 2)) for c in X]#パワースペクトルへ変換[dBSPL]  # 振幅スペクトル
        phaseSpectrum = [np.arctan2(int(c.imag), int(c.real)) for c in X]    # 位相スペクトル

        #----ピーク取得-----+
        #ピーク値をストックする
        lvl, frq = get_peak(amplitudeSpectrum)
        for x in xrange(len(lvl)):
            all_lvl.append( lvl[x] )
            all_frq.append( frq[x] )
            #print x

        #------図表書き出し------#
        #all_lvl,all_frqにピーク値をストック
        #最終ループですべて書き出す
        if file == len(files)-1 : #x ==( len( lvl)-1 ):
            plt.figure(figsize =(4,3))
            #subplot(211)
            plt.scatter(all_frq, all_lvl,s = 2,c = 'b', marker = 'o')
            axis([0, fs/2, 0, 100])
            xlim([10,500])
            ylim([0,50])
            xscale("log")
            plt.savefig('test.png'     ,format = 'png', dpi=300)
            #plt.savefig('peakpoint.png', format = 'png',dpi = '300' )
        
        del lvl[:],frq[:] #値をクリアして次のループへ

ディレクトリに沢山の時系列データ(今回は**.csv)がある場所で使うと,個数分だけループをはじめて,1つ読み込んでFFT処理,ピーク値をプロット,また1つ読み込んでFFTしてピークとってきてさっきの絵にプロット…というプログラム

(実際の挙動は一回一回プロットしてるとアホらしいのでプロットすべき点をストックさせて一気にプロットする)


何が便利かというと,ある室に関して,沢山の受音点データがある時これで周波数特性のピーク値をプロットするとピークがどの周波数にあるか点で見ることができる.

f:id:coupe-glass:20170514000526p:plain:w400:h300
今回は361点分の受音データ(全て同じ室)を処理したんだけどこんな感じで30Hz付近にピークが集まってる=定在波が立つかなっていうおおよその検討がつけられるようになった.


半年ぐらい前に作るべきシステムだったんだけどPythonの知識がなかった….