科捜研の人

女房が実家から大きな包みを持ってきた。いぶかしげに何それ？ って、聞いたら、DVDも見れる 14インチサイズのテレビでしたよ。

なんでもそのTV、実家に貰われてきたそうなんだけど、番組が写らないんで、何が悪いかおいらに 調べてみて欲しいとの事。

5cm X 15cmぐらいなモバイルアンテナが付属してた。早速それを取り付けて、この地に 馴染むようにスキャンさせたよ。でも、何も写らん。アンテナを高い所に出して、外にあるさかなの 骨アンテナと同じ向きにしたけど、写らなかった。

しょうがないので、壁に来てるTVコンセントに繋いでスキャンさせたら、やっとこの地に馴染んで 全チャンネルが受信出来るようになった。

早速実家へ持って行って、電波の入りやすそうな所で試したら、１局受信出来た。後はARDF宜しく、 電波が来てそうな方向目指して、右往左往しましたよ。それでやっと全チャンネル受信出来る ようになった。やれやれ、ARDFなんて初めての体験だったわい。

ちなみに携帯アンテナは、DUA-300 だった。ゲインがほぼ0って事は、中にプリント基板が一枚入ってて、アンテナパターンが書いて あるだけだな。 この地の屋外アンテナは20エレメントぐらいで13dBぐらいある事を思えば、あんな 板アンテナでよく写ったな。

音量

TV大好き人間の女房は炬燵に入って、リモコンでZappingしてる。CMが入る たびにやってるよ。どこのご家庭でもそうなんですかね？ おいらはアコギなCMが入るんで、 ほとんどTVを見なくなっちゃったぞ。

CMになると音量が大きくなるのもうざいな。そのたびにリモコンの音量ボタンを調整してるよ。 このCMの音量増大は世の不評を買ってるみたいで、今年の10月からは改善されるとか。そんなの 明日から出来るじゃん。やらないのははTV局がスポンサー様に気兼ねしてんのね。

CMでっせって信号で録画を中止するレコーダーをどこかのメーカーが作っていたはずだけど、 これ業界から圧力がかかって発売中止になっちゃったんだよな。

CM時の音量って、CM以外の時とどのぐらい音量が違うんだろう？ 騒音計で測れるのかな？ CMなんて騒音以外の何物でもないから、きっと測れるのだろうな？

音量って時々刻々変化してくから、一定時間毎の音量を求める事になるのかな？パソコンで 録音して処理出来そうだな。

音声信号って交流だから、半波整流なり全波整流して、それを積分すればいいのか。積分時間が 計測時間になるんだな。

半波整流は、正の信号だけ取り出せばいいし、全波整流なら、絶対値を取ればいいんか。積分は 積算に置き換えれば済むか。

まてよ、耳には音の周波数による感度差があるから、それじゃ聴覚特性を無視した計測になっちゃうな。 どうすればいいのかな？

FFTして周波数毎に積算かな。結果表示は？ ふと思ったんだけど、 科捜研の人が、犯人からかかってくる電話を声紋分析してたな。あれは縦軸が周波数で横軸が 時間だったかな。音の強弱が色の違いになって表示されたはず。そんなソフトって簡単に 手に入る？

先生に聞いてみよっと。

Pythonによる音声信号処理

Pythonで音声信号処理 で サウンドスペクトラムってのが ぴったりだ。

このページで紹介されてるのは、録音済みのwavファイルからスペクトラムを 造ってくれるやつ。簡単に試せて面白い。

import wave
from pylab import *

if __name__ == "__main__":
    wf = wave.open("aaa.wav", "rb")
    data = wf.readframes(wf.getnframes())
    data = frombuffer(data, dtype="int16")
    length = float(wf.getnframes()) / wf.getframerate() 
    
    N = 512
    
    hammingWindow = np.hamming(N)

    pxx, freqs, bins, im = specgram(data, NFFT=N, Fs=wf.getframerate(), noverlap=0, window=hammingWindow)
    axis([0, length, 0, wf.getframerate() / 2])
    xlabel("time [second]")
    ylabel("frequency [Hz]")

    show()

ついでなので、specgram がどうなってるか見てみる。ソースは、 matplotlib/mlab.py に有る。冒頭部分にコメントが入ってた。

def specgram(x, NFFT=256, Fs=2, detrend=detrend_none, window=window_hanning,
        noverlap=128, pad_to=None, sides='default', scale_by_freq=None):
    """
    Compute a spectrogram of data in *x*.  Data are split into *NFFT*
    length segements and the PSD of each section is computed.  The
    windowing function *window* is applied to each segment, and the
    amount of overlap of each segment is specified with *noverlap*.

    If *x* is real (i.e. non-complex) only the spectrum of the positive
    frequencie is returned.  If *x* is complex then the complete
    spectrum is returned.

    %(PSD)s

    Returns a tuple (*Pxx*, *freqs*, *t*):

         - *Pxx*: 2-D array, columns are the periodograms of
           successive segments

         - *freqs*: 1-D array of frequencies corresponding to the rows
           in Pxx

         - *t*: 1-D array of times corresponding to midpoints of
           segments.

詳細は、下請けメソッドの _spectral_helper に任せているのね。詳細を見てくと、NFFT毎に 区切ってfftしてるな。それが終わると、X軸（時間方向の目盛り）とY軸（周波数方向の目盛り） を計算してるんか。

折角なので、ちとどんな風になってるか調べてみる。まずは対象となるaaa.wavの素性をsox （次の節に説明あり）で調べておく。

C:\homes\CLJ>sox --info aaa.wav

Input File     : 'aaa.wav'
Channels       : 1
Sample Rate    : 8000
Precision      : 16-bit
Duration       : 00:00:01.53 = 12260 samples ~ 114.938 CDDA sectors
File Size      : 24.6k
Bit Rate       : 128k
Sample Encoding: 16-bit Signed Integer PCM

統計情報に、大雑把は周波数が出て来るなんて、面白いな。

C:\homes\CLJ>sox aaa.wav -n stat
Samples read:             12260
Length (seconds):      1.532500
Scaled by:         2147483647.0
Maximum amplitude:     0.339142
Minimum amplitude:    -0.371918
Midline amplitude:    -0.016388
Mean    norm:          0.034040
Mean    amplitude:    -0.000000
RMS     amplitude:     0.055354
Maximum delta:         0.528748
Minimum delta:         0.000000
Mean    delta:         0.034994
RMS     delta:         0.057137
Rough   frequency:         1314
Volume adjustment:        2.689

上記のスクリプトでspecgramを抜けた所で、返値を調べてみる。

(Pdb) bins
array([ 0.032,  0.096,  0.16 ,  0.224,  0.288,  0.352,  0.416,  0.48 ,
        0.544,  0.608,  0.672,  0.736,  0.8  ,  0.864,  0.928,  0.992,
        1.056,  1.12 ,  1.184,  1.248,  1.312,  1.376,  1.44 ])
(Pdb) freqs
array([    0.   ,    15.625,    31.25 ,    46.875,    62.5  ,    78.125,
          93.75 ,   109.375,   125.   ,   140.625,   156.25 ,   171.875,
            :
        3843.75 ,  3859.375,  3875.   ,  3890.625,  3906.25 ,  3921.875,
        3937.5  ,  3953.125,  3968.75 ,  3984.375,  4000.   ])
(Pdb) size(bins)
23
(Pdb) size(freqs)
257

binsってのがX軸でfreqsはY軸の目盛りか。データは8KHzサンプリングなんで、FFT すると、4000Hz以上は折り返しで不用になるんだな。

(Pdb) 12260 / 8000.0
1.5325
(Pdb) 12260 / 512.0
23.9453125

X軸の単位は秒だけど、データ数から算出。FFTの回数も１回のデータ数から算出してんだな。

(Pdb) pxx
array([[  1.96987443e-02,   1.05299310e-01,   1.21728279e-01, ...,
          5.29375972e-02,   3.12171009e-03,   5.46604318e-03],
       [  2.40170120e-02,   4.58577788e-02,   1.16091109e-01, ...,
          2.91492652e-02,   2.22824779e-02,   4.67810856e-02],
       [  2.61483267e-03,   1.02024403e-01,   1.85373853e-01, ...,
          4.12349123e-02,   1.72837381e-02,   1.27102490e-01],
       ...,
       [  3.18807711e+00,   8.69322316e-01,   4.21885091e-01, ...,
          4.95199933e-01,   6.98053060e+00,   3.14729810e-02],
       [  5.11484378e-01,   1.67718000e+00,   9.63667184e-01, ...,
          3.47356476e-01,   1.38390202e+01,   1.47733175e-01],
       [  2.96100749e-01,   2.75793445e-01,   3.54159019e-01, ...,
          1.24707272e-01,   9.89525629e+00,   8.61189585e-02]])
(Pdb) size(pxx)
5911
(Pdb) size(pxx[0])
23
(Pdb) 5911 / 23
257

pxxは２次元配列か。一番低い周波数から高い周波数方向に配列が257個連なっている。そして 個々の配列内は、時間毎のスペクトラムのピークが並んでいる。

後は、show()で、データの大きさによって色に置き換えて表示してるんだな。納得しましただ。

後、こちらのページ も面白いぞ。

sox

上のページで紹介されたsoxってのが、なかなか綺麗な図を 書いてくれるな。

sox aiueo.wav -n spectrogram

これで簡単にpngファイルが出来上がる。soxって本来こんな使い方をするものじゃないのね。 SoXの基本的な説明とか SoX - Sound eXchange に関するメモに 書いてあるのが正統な使用法なんで、お間違えなき様。

これだけじゃ、soxにニアミスしただけでつまんないからソースを眺めてみる。もろにspectrogram.c が有るんだけどちょいと複雑。これはもう得意のgdbの出番ですな。普通にconfigureしただけじゃ spectrogramが使えないので、

CPPFLAGS=-I/usr/local/include LDFLAGS=-L/usr/local/lib ./configure   --enable-debug

した。で、以前に造った vim + Pyclewm を使ってsoxを呼び出そうとしたら、srcの直下にある soxは、シェルスクリプトだったよ。何が起こっているか、 sh -x ./sox したら、本体が、 src/.libs/sox内に隠されていたよ。従って、本体を直接呼び出してしまった。

soxって、カメレオンみたいに名前を変えて呼び出せるんだな。playとかrecとかsoxiとか。だから シェルスクリプトで包んで、分かり易くしてんのか。面白いな。

spectrogram.cのstartにBPを置いてgdbを走らせてみた。

  Breakpoint 2, main (argc=4, argv=0xbfbfe6a8) at sox.c:2688
  (gdb) continue
  [Inferior 1 (process 4011) exited normally]
  (gdb)

けど、breakしない。どうやらお目当てなやつは、サブプロセスとして実行されちゃっててgdbの 制御範囲外になってるようだ。こういうのはどうやって追跡したらいいのだろう？ 後で調べて みるかな。

audacity

上記のsoxのGUI版が、audacityだ。 FreeBSDでもpkgが有るので入れるのはやぶさかでなないけど、今回は見送った。Win用が入って いれば事が足りるからね。

録音を始めると、トラックに音の波形が出てきたり、レベルメーター（バーだけど）がほいほいと 動いたりして、ああ録音してるんだなって気にさせられる。

音声トラックの所に付いている逆三角形のアイコンをクリックすると、スペクトラム表示に 切り替えられるんだけど、サンプリング条件によっては、表示が大幅に遅延したり、全く 表示されなかったりで、ちとこつがいるようだ。

その他の機能として、波形を切り取ったり、エフェクトをかけたり、いろいろな形式で録音 データをセーブ出来たりする。至れり尽くせりなアプリですな。HDDに入れといて損は無いと 思うよ。

SFS/WASP

audacityもいいんだけど、スペクトラムに注目するんなら、やはりそれ用に造られたものの 方が操作も単純でいいと思うよ。

探してみると、売り物のやつとか（その試用版は）沢山見つかるんだけど、中々フリーの ものって無いのね。で、やっと見つけたのが、 Windows Tool for Speech Analysis ってやつ。英国の大学って事は産学協同って事で、スコットランドヤードとタイアップして るのかな。Human Voices and the Wah Pedal こんな研究も進んでいるしね。

録音したやつをスペクトラム表示したり、フォルマント表示したりで、より捜査用って感じが するな。きっと科捜研の人も使ってて、方言のわずかな違いを認識出来たりして。おー怖。

SFS/RTGRAM

上のやつは、録音してから解析って事なので、ちとリアルタイム性に欠ける。本当にマイクから 入った信号をその場でスペクトラム表示出来るのは無いかと探したよ。そしたら、 Windows Tool for Real-time Speech Spectrogram Display ってのが見つかった。上のやつと姉妹編なのね。

折角だから、らじる★らじるを信号源にして、リアルタイム に表示してみた。これでニュースでも表示しながら、声紋から担当アナウンサーを推測するって のが面白いかも知れんな。嗚呼、アナウンサーを特定するなら、 ラジコの方がいいかな。

ハム用のトランシーバー FT-817 をつないで40mを聞いてみたら、CWの線が３本も４本も表示 されて、目でも分かる混雑ぶりでしたよ。これ、横に流れるウォーター・フォールだな。 正に貧乏人御用達であります。

おまけ

Pythonのコードをさっと見る時、余計な物が目に入ってしまう。余計な物が折り畳まれていると 便利だ。そんな時は、Fold python code nicely and toggle with one keystroke を入れておくと良い。

ファイルを読み込むと、メソッドが折り畳まれた状態で表示される。見たい所へ移動して、f キー を押すと展開される。もう一度 f キーを押すと畳まれる。fの代わりにFキーを使うと、ファイル内の 全てのメソッドを畳んだり展開したり出来る。これ、emacsにも有ったな。

もう一つあって、vimの中からpydocを引けるやつ。 Python documentation view- and search-tool (uses pydoc) から落としてくれば良い。:Pydoc word とかやってもいいし。引きたいwordの所へカーソルを 合せてから、<leader>pw としても良い。