鳥の声の分析はオープンソフトの"R"や"Amadeus Pro(MacOSX)"を用いて来た。
$ vi BPF.m
それを3Dで見てみるとこんな感じでうまく行っている様です。
これはこんな自前の短時間FFTで表示しました。もちろんこちらの本を参考にして。
{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
所がである。声紋（ソナグラム）をもっと詳細に分析したいと思ったのだが、なかなが歯が立たない。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番目は高い方のカットオフ周波数。
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
コメント