スキップしてメイン コンテンツに移動

Mathematicaで鳥の声解析　(BPF:Band Pass Filter)

鳥の声の分析はオープンソフトの"R"や"Amadeus Pro(MacOSX)"を用いて来た。

所がである。声紋（ソナグラム）をもっと詳細に分析したいと思ったのだが、なかなが歯が立たない。Amadeusはソナグラムはもちろん表示できるが、瞬時周波数の表示が面倒なのだ。Audacityは使い勝手が悪いし。

そこで買ったままになっていたMathematica 9.0を再び起動させることにした。今はver.11にアップデートしませんかと催促のメールが来るが今のところ放置。

という訳で、MacOSXのTerminalからMathematicaで作ったプログラムをフィルタのように使いましょうというお話。

コマンドラインからMathematicaを起動してパラメータも渡し、結果をUNIXコマンドみたいにパイプに渡すことができる。

今回作ったのは、バンドパスフィルタ（BPF）のフィルタ。

MacOSXのTerminalからこんな感じで使います。

$ ./BPF.m input.wav output.wav 1000 14000
1番目のパラメータは自分のコマンド（BPF.m）だとして、2番目は入力音声データ、3番目は出力データ、4番目は低い方のカットオフ周波数、5番目は高い方のカットオフ周波数。



$ vi BPF.m
   1 #!/Applications/Mathematica.app/Contents/MacOS/MathematicaScript -script
  2 
  3 (* BPF.m : BandPassFilter *)
  4 If[ Length[$ScriptCommandLine ] < 4 ,
  5  Print["programname=BPF.m (mathematica script ver.1.3):入力波形（モノラル）にバンドパスフィルタをフィルタを掛けて波形出力する"]
  6  Print[" Usage:Ver1.3 Fourier変換のための個数を偶数化する"]
  7  Print[" Usage: カットオフ低周波数(fcl)以下とカットオフ高周波数(fch)以上のスペクトルを0にして波形に戻す"]
  8  Print[" Usage: 出力はDesktopになされる"]
  9  Print[" Usage: BPF.m freq_low freq_high inputfile"]
 10  Print[" Usage: BPF.m 6000 14000 a150428_065_10sec.wav"]
 11  (*  Print[" Usage: BPF.m freq_low freq_high inputfile outputfile"] *)
 12  (* Print[" Usage: BPF.m 6000 14000 a150428_065_10sec.wav ~/Desktop/a150428_065_10sec_BPF.wav"]*)
 13 
 14  Exit[]
 15 ]
 16 
 17 cutLowFreq      = ToExpression[$ScriptCommandLine[[2]]]//N;
 18 cutHighFreq     = ToExpression[$ScriptCommandLine[[3]]]//N;
 19 inputfile       = $ScriptCommandLine[[4]];
 20 Print["inputfile=",inputfile]
 21 (* outputfile   = $ScriptCommandLine[[5]]; *)
 22 outputfile      = FileNameJoin[{"~/Desktop", inputfile}]
 23 ops             = Import[inputfile, "Options"]; (* 各変数が取り出せる *)
 24 wave            = Import[inputfile, "Data"]; (* 音声データの読み込み *)
 25 sampleRate      = Replace["SampleRate", ops]; (* Print["sampleRate=",sampleRate]; *)(* サンプリングレートの設定 *)
 26 
 27 myBPF[data_List, fcl_, fch_, samplingRate_: 48000.] :=
 28   Module[{wave,ml,fw, fmin, bpfFilter},
 29     ml=2*Round[Length[data]/2];
 30     wave=data[[1;;ml]];
 31     fw = Fourier[wave];
 32     fmin = samplingRate/ml;
 33     bpfFilter[f_] :=  (* バンドパスの周波数設定 *)
 34       Which[
 35         f < fcl , 0,  (* fcl未満は0に設定 *)
 36         f > fch , 0,  (* fch超は0に設定   *)
 37         True    , 1   (* その他は1に設定  *)
 38       ] ;
 39   bpf1 = bpfFilter[#] & /@ Table[i*fmin, {i, 1, ml/2}]; (* Fourierは半分で折り返すので下半分を用意 *)
 40   bpf  = Join[bpf1, Reverse[bpf1]]; (* 上半分をリバースでJoinする *)
 41   InverseFourier[fw bpf] (* 入力スペクトル(fw)にBPF(bpf)を掛けて、それを逆Fourier *)
 42   ]
 43 
 44 waveBPF = Re[myBPF[wave, cutLowFreq, cutHighFreq, sampleRate]]; (* モノラルのみ対応 *)
 45 snd = Sound[SampledSoundList[waveBPF, sampleRate]]
 46 Export[outputfile, snd]

もちろん行頭の数字はソースには入っていません。また、実行プログラムにするため、
$chmod +x HPF.m
は忘れずに。1行目がお呪い。

このフィルタは指定されたカットオフの低周波数(fcl)以下と高周波数(fch)超をばっさり”0”にします。それ以上は”1”でステップ関数。これを音声スペクトルに掛け合わせしています（ステップ型です）。このためサイドローブが出る事も心配ですが、欲しい周波数と離れていればまあ大丈夫でしょう。

 それを3Dで見てみるとこんな感じでうまく行っている様です。

これはこんな自前の短時間FFTで表示しました。もちろんこちらの本を参考にして。

waveR = waveBPF[[1]];
ffts = Table[Fourier[waveR[[t ;; t + 1024 - 1]]] , 
  {t, 1, Length[waveR] - 1024, 1024/2}];
lffts = Abs[ffts[[All, 1 ;; 512]]]; 
ListPlot3D[lffts, Mesh ->; None, 
 PlotRange -> All,
 ImageSize -> 500,
 SphericalRegion -> True,
 ViewPoint -> {-1, 1, 0.75},
 Ticks -> {Function[{min, max},
 Table[{i, i*30}, {i, 0, Floor[max], 100}]],
   Function[{min, max},
   Table[{i, i*0.33 // Round},
 {i, 0, Floor[max], 10}]], Automatic} 
]

参考図書：




Ver 1.3：音声データの個数が奇数だとfilterとFourier変換の個数が違うのだろう。偶数にする。
Ver 1.2 ：やっぱりステレオはうまく行かないのでモノラルにしています
Ver 1.1 : ソースはこんな感じ。入力ファイルはモノラルでもステレオでもOK

ラベル:

コメント

コメントを投稿

このブログの人気の投稿

瞬間日記のデータをday oneに移動させるのだ

瞬間日記をiPodを購入した２０１０年から使ってきたが、day oneが人気なので使ってみた。なかなかよかった。コマンドライン(CLI:Command Line Interface)から入力できるし、クラウドで同期できるのもいい。見た目もきれい。安定しているのもいい。 という訳で、瞬間日記からday oneにデータを引っ越しすることにした。 （このページで半角の ＆ｇｔ； と表示されたら半角の”＞”記号だと思ってください。HTMLの仕様のインプリバグでしょう。） 特定の日記アプリ間のデータの移動とは人生で一回しかないだろうから汎用性など考えずに間違えないように慎重にすることが大切だ。確かめて確かめて。このページはそんな備忘録である。また、瞬間日記のデータをday oneに移動させたという記事も見つけられなかったので誰かのお役に立てるかもしれない。でも自己責任でね。 瞬間日記は独自フォーマット(.dat3)や標準的な(.sqlite)、(.csv)でもMacにバックアップできる。使い勝手のいいcsvファイルは本文だけで写真情報が欠落している。瞬間日記のデータを引き上げるのは色々試したが、テキストばかりでなく写真も移動させたいの結局メールでバックアップすることにした。瞬間日記側で自分にひと月毎のバックアップをメールする。３年半分。ソチオリンピック開会式を見ながら作業でもれなく終了。 これは瞬間日記が３０枚しか１つのメールに写真を貼付できないので、ひと月分だけせっせとメールする。 届いたメールはOSX側のmail.appで開くが、~/Libraryに保存されているelmsでなく、単にテキストファイルでmail.app側の機能を使って出力することにする。 瞬間日記のデータを取り出し、day oneにエントリできるくらいの粒度のファイルにするまでの流れ： iPodの瞬間日記のひと月分に区切ったバックアップデータを自分のメールアドレスに送信する OSXのmail.appで受信したメール全部（件名"MomentDialy"で始まるメールを一つのメールボックスに束ねる。例えばMomentDiaryというメールボックスに全部入れる mail.app上のメールボックスMomentDiary内の全部のメールを選択する。つまり⌘+a。 ...
コメントを投稿
続きを読む

鳴き声からの識別(How to identify the bird from its sound. )

鳥はその暮らしの中で様々な音を出します。 朗々と鳴く囀り、冬の地鳴き、飛ぶ時の羽音、何かを突っつく時の音、翼を体にぶつけて出す母衣（ほろ）うち、虫を探しているのか枯葉をひっくり返すときにでるカサコソと言う音。その中にあって、自己を主張しコミュニケーションとしての囀りや地鳴きの声は鳥の名前を探り当てる大いなる情報です。 鳥の鳴き声を聞いて、なんて言う鳥が鳴いているんだろう、とバードウォッチングの最中に思う事はしょっちゅうです。ベテランのバードウォッチャーはこの声は○○と教えてくれます。では、聞く相手も居ないとき、バードウォッチングしている仲間の誰も分からないとき、鳥の声からどう鳥の種名を探り当てる事ができるでしょうか。このページではそんな事を考察したいと思います。 参考にするレファレンスCDはこの4冊：蒲谷1996（蒲谷2004）、上田1998、松田2016、植田2014、このページの下にリンクを張りました。この他にも、 バードリサーチ鳴き声図鑑 もあります。 xeno-canto も有名です。 それでは、鳥の鳴き声から種名同定までの幾つかの方法を見て行きましょう。 初級編 方法1）探鳥会で人に聞く 当たり前ですが、知っている人に聞くのが一番手っ取り早いです。 バードウォッチングの初心者は見るモノ聞くモノすべて新鮮で教えてくれる人もたくさんいます。探鳥会に参加してその場所で年中見られる普通種をまず覚える事になります。これが”ものさし”となってこれと同じか違うかが分かれば、聞き分けられる数が多くなってきます。 方法2）さえずりナビを試す。 バードリサーチが”さえずりナビ”をスマートフォン向けに作っています。 家の周りで，ハイキングに行ったとき，鳴いている鳥の種名を知りたくなることがあります．「さえずりナビ」は，場所・季節・声のタイプなどから，その鳴き声の鳥を検索してくれます．姿や形の専門知識が無くても検索ができ，その鳴き声を聞くことができます． あまり使った事無いのですが、twitterではまあまあのコメントが多いです。 中級編 方法3）録音して似ている声をCD/WEBで調べる 録音できれば後で調べる事ができます。ICレコーダも小さく安くなってきました。持っていなくても携帯電話やコンパクトカメラの動画機能でとりあ...
コメントを投稿
続きを読む

電解コンデンサの容量抜けを測ってみた

中古オーディオ屋でスーパーウーハーを買ってから7〜8年経つだろうか。 1989年に発売されたONKYOの SL-10 という機種で、現役で商品が店頭に並んでいた時にピアノ曲に深みが増すことにちょっとした感動を覚えて中古屋で見つけた時2万円で即購入したのだ。 スピーカのエッジのクズがエンクロージャのダクトから出てきて開けてみたら見事エッジ全体が下手っていた。これを奇麗にして、ウレタンエッジに交換することにした。写真1、と写真2がそのBefore/Afterだ。 写真1　SL-10のスピーカのエッジを交換してみた （エッジを取り除いた状態） 写真2　SL-10のエッジを交換したスピーカ また、SL-10が勝手に自己発振する様になっていたのでコンデンサの容量が抜けたと想像してこれらを交換することにした。これは電解コンデンサの容量が抜ける事でフィードバック回路の時定数が変わってしまいネガティブフィードバックがポジティブフィードバックに位相が回って知ったのでは予想した。 電解コンデンサを交換したSL-10のコンデンサ 交換した電解コンデンサ達（容量は抜けていなかった）  マルツ電波で発注しておいたオーディオ用のコンデンサを付け替えた。ただ、それだけでは詰まらないので交換したコンデンサを測ってみた。使った測定器は卓上テスタ（GBW 9000A）の付属機能。最大20μFまでしか測れなかったが、結論からすると交換した電解コンデンサの容量抜けは無かった。 測定例 定格（μF） 　　実測（μF） 10　　　　　　　10.19 10　　　　　　　10.00 4.7　　　　　　　4.78 3.3　　　　　　　3.27 1.0　　　　　　　1.01 0.47　　　　　　　0.46 というわけで直しても自己発振は直らないかも。明日組立ててみよう。
2 件のコメント
続きを読む