音響カメラが欲しいという話

皆さんこんにちは。鳥さんのおちりを追いかけている者です。最近暑くなってきて長時間外に出るのもしんどくなってきているのではないでしょうか？
少なくとも私はそうで、そうなると鳥さんを探す時間も減ってきてしまいます。
さらに野鳥の撮影をされない方からは賛同を得られないかもしれませんが、夏が近づき木が葉でもっこもこになっているとどこに鳥さんがいるのかわからないのですね。

という事で鳥さんの声が聞こえたら速攻でそちらの方向にカメラを向けないといけないわけです。
ただし初心者からすると鳴き声はするけどどこにいるかわからん…
となってしまいますので、いわゆる音響カメラが欲しいなーと思ってしまったわけです。

しかしながら音響カメラ自体があまり有名な業界ではないのでお手軽に買うことができない…という事でどうにかして作ることができないかと考えたわけでございます。
当然思い付きでは出来なさそうですので、とりあえず自分で考えたシステムで何とかなりそうかというところまで考えてみました。

１，理想的な条件とは

普通音響カメラには複数個のマイクが取り付けられており、それぞれのマイクに到達した音波のずれ量からどこから音が出てきているかを計算することで位置を特定します。このため当然ながら音の波形はどんなにずらしても完全に一致することはありません。
しかしながらまずは同一の音声データを用意し、それらを数ミリ秒ずらし２つのデータにします。元のデータをf(t)としたときにf(t+dt)をプログラムに入れてdtを求めることができるかという事をまずはやってみようという事にしました。
これが理想的な条件というわけですね。ノイズもないし答え合わせも簡単なはず…

２，どーやって時間差を検出するのかという問題

とりあえずPCとかのパワーに任せて解析をするとしたらこんな風にやってみましょう
A.とりあえずFFTを実施
B.特定の周波数帯のピークを計算する
C.Bで取得したピークの周波数fの位相を計算する
D.２つのデータの位相差Φを計算しdt=Φ/f
完璧なプランですね。

３，という事で実際にやってみたところ

ソフトの実装については省略して結果だけを示します。STFTのように全音声データを細切れにして上の手法で時間差を解析してみました。利用した音声データとかはこんな感じ。

音声データのサンプリングレート：48kHz
FFTのサンプル数：2048
ずらしたフレーム数：6=0.000125秒

という事でこのデータを解析し、頻度分析をしたところこんな感じの結果になりました。ピークをとっている位置が無事に0.00012くらいになっておりました。

理想的なデータならそこそこの精度で計算できているか？

今後の予定

とりあえず次は実際にマイク２つを用意してデータを取ってみたいと思います。うまくいかなかったら更新はありません。

ではまた

import pydub
import numpy as np
from pydub import AudioSegment
import matplotlib.pyplot as plt
import os
import matplotlib.pyplot as plt

if __name__ == "__main__":

    song = AudioSegment.from_wav(r"C:")
    sample_rate=song.frame_rate
    print("サンプルレートは"+str(sample_rate)+"Hzです")
    
    song_array = np.array(song.get_array_of_samples())
    song_ch1 = song_array[::2]
    song_ch2 = song_array[1::2]
    
    data_size=song_ch1.size
    wsize = 2048
    delta=6
    f_range_low=400
    f_range_hi=2000
    
    #周波数の決定
    freq=np.fft.fftfreq(wsize,1/sample_rate)   
    idx_lo=int(f_range_low/(sample_rate/wsize))
    idx_hi=int(f_range_hi/(sample_rate/wsize))
    freq=freq[idx_lo:idx_hi]
    freq_rad=freq*2*np.pi
    
    ret=[]
 
    for i in range(0,data_size-wsize,200):
        data=song_ch1[i:i+wsize-1]
        data_delta=song_ch1[i+delta:i+wsize-1+delta]
        
        #FFTを計算し、特定の周波数帯のピークを探す
        
        fft_result = np.fft.fft(data)
        amp = 2*np.abs(fft_result / wsize)
        deg = np.angle( fft_result )
        
        fft_result_delta = np.fft.fft(data_delta)
        amp_delta = 2*np.abs(fft_result_delta / wsize)
        deg_delta = np.angle( fft_result_delta )
        
        #解析に利用するデータの切り取り
        amp_trim=amp[idx_lo:idx_hi]
        deg_trim=deg[idx_lo:idx_hi]
        amp_delta_trim=amp_delta[idx_lo:idx_hi]
        deg_delta_trim=deg_delta[idx_lo:idx_hi]
        #ピークの位置計算
        idx_pk=np.argmax(amp_trim)
        idx_delta_pk=np.argmax(amp_delta_trim)
        dphi=deg_trim[idx_pk]-deg_delta_trim[idx_delta_pk]
        dt=dphi/freq_rad[idx_pk]
        #時間の誤差
        ret.append(dt)
        

n, bins, patches = plt.hist(ret, bins=100, range=(-0.001, 0))  

# ヒストグラムの表示  
plt.show()

print("end")