spectrogram (2)
前回は思わぬ事から、ラズパイの石の裏の顔を知った。裏の顔なんで、まだ余り世間一般には知られていないようだ。
複雑過ぎる石は、体に毒なんで優しい石を調べてみるか。
// 取得した値を電圧に変換 (0〜5000mV) tempValue = map(tempValue, 0, 1023, 0, 5000);
mapの鏡みたいな関数だなあ。どなたにもご理解頂けます。所で、アルディーノのウノが有名みたいだけど、石は何処のやつ？
以前にcw decoderを移植した時、millsって関数があった。50日でオーバーフローするから、気を付けてねってやつ。何Bitのカウンターを実装してるか予想してみる。
octave:3> log2(1000 * 3600 * 24 * 50) ans = 32.008
32Bitのバイナリーカウンターだってさ。
Arduino は小さなコンピュータ・システムを経由して、いきなり石の仕様書に辿りついた。
ざっと見、32Bitのカウンターなんて実装されていないぞ。でも、内蔵のA/Dコンバーターを使って、チップの温度を表示出来る事を知った。面白いな。
mic
今まで声紋の代表として、カラオケCDから取り出した物を使っていた。 けど、この後に及んで、 自分の声紋を見たくなった。 声紋ってものの本によると、1962年にベル研のKerstaが、spectragramによる話者認識の可能性を発表した事に端を発しているそうだ。
そんな訳で、録音機が欲しい。目の前にあるじゃん。Windows10のボイスレコーダーで、録音。 取り合えず、短いフレーズで「開けごま」。そこまではいいんだけど、結果はWindows仕様の音形式。まったく、田舎仕様なんで困っちゃうなー。何とかならないかね。
(base) cent:~$ file me.m4a me.m4a: ISO Media, MPEG v4 system, version 2 (base) cent:~$ soxi me.m4a soxi FAIL formats: no handler for file extension `m4a'
CentOSは知っていたとな。でも、soxさんは、そんな物扱いに困りますってね。そこで例のやつ登場。
(base) cent:~$ ffmpeg -i me.m4a me.wav
:
Stream #0:0 -> #0:0 (aac (native) -> pcm_s16le (native))
Output #0, wav, to 'me.wav':
Metadata:
major_brand : mp42
minor_version : 0
compatible_brands: mp41isom
album_artist : ボイス レコーダー
ICRD : 2019
生意気に、ステレオ録音だったよ。wavに変換すると10倍に膨れ上がった。それから48kサンプリングと言う豪華さ。こういう音じゃないと、コルチナは聞き分け出来ないんだな。
でも、パソコンに向かって喋るの性に合わないな。なんだか、コルチナだかシリだかアレックスだか、外人に指示してるみたいで。日本人なら、おーいお茶って叫びたい。
以前使ってたマイクがあったな。確か胸に付けるピンマイクって言うやつ。マイクを口の前に置くと、吹くと言う現象に見舞われるそうです。パピプペポ系の発声は、マイクに息が突風のごとく押し寄せ、思わぬ雑音となるそうです。プロはその防止の為に、マイクの前に暴風幕を設置するんだそうです。
専用装置が有るそうですが、無ければ簡単に自作出来るとか。ストッキングを刺しゅうの時に使う枠で固定すれば良いそうです。
ヘッドフォンとマイクがセットになった物なら、マイクを口角の脇にセッティングするのが良いそうです。ピンマイクなら、胸に付けるだけなんで、そんな面倒は不要。
rec myvoice.wav trim 0 10
とかすれば、10秒間だけ録音できます。おっと、Windows10のパソコンは、スカイプ用の仕様みたいで、ジャックがマイクとヘッドフォンの一体型になってた。つまんないな。
AudaCity
そんな訳で、debian機の方で確認。どうせなら、DISKの肥やしになっているaudacityに登場願おう。
録音ボタンを押したら、トラックが２本出てきて、録音が始まった。ボイスを確認しながら録音するって、一端のプロになったみたいで気持ちいい。
どんな機能が有るか一応確認しとくか。
マルチトラック編集に対応した、高機能オーディオエディタ！「Audacity」
praat
これ、たまたま見つけたもの。フォルマントとかを認識できるみたい。 debianにも有るので入れて、起動だけは確認した。下記を参考にやっていけばいいのかな。
at golang
golang方面のaudio関係はどうなってるかと思って調べてみた。有志が提供してたね。 声紋を出すなら取り合えずFFTが有ればなんとかなる。
有ると分かれば安心して、次に行けるな。
at octave
と言う事で、懐かしいoctaveです。貧者のMATLABです。
Scientific Programming Language GNU Octave
高専の教授が入門資料を公開されてました。
信号処理分野のためのGNUOctave/MATLAB/Scilab入門
そして、サンプルも提供されてました。ありがたい事です。
Demo Programs of Signal Processing
on debian
debianに古いoctave(4.0.3)が入っていたので、スペクトラグラムしようとした。声紋は普通の人が使う物ではないので、頑張っていれる。
octave:6> pkg install -forge signal octave:7> pkg load signal
これで、HOME直下にoctaveが出来て、その中にパッケージが入った。使う時はloadをお忘れなく。デモ出来るみたい。
octave:13> demo specgram 1 specgram example 1: Fs=1000; x = chirp([0:1/Fs:2],0,2,500); # freq. sweep from 0-500 over 2 sec. step=ceil(20*Fs/1000); # one spectral slice every 20 ms window=ceil(100*Fs/1000); # 100 ms data window ## test of automatic plot [S, f, t] = specgram(x); specgram(x, 2^nextpow2(window), Fs, window, window-step);
白けたグラフしか出てこなかった。グラフのメニューをつつくと、声紋が出てきた。どういうこっちゃ？
on OpenBSD
じゃ、ウブとかFreeBSDって手が有るけど、java,perl,pythonとか、色々余計な物が付いてくるので敬遠。OpenBSDに4.4.1が有ったので入れた。余り余計な物はついてこない。
少し資料を漁ってみた。
お手製の声紋作成スクリプトがあったので、頂いた。
ob$ cat spectrogram.m
function [S,time,frequency]=spectrogram(s,fs,bits,window_size,shift_size)
length_of_s=length(s);
number_of_frame=floor((length_of_s-(window_size-shift_size))/shift_size);
N=window_size;
x=zeros(1,N);
w=hanning(N);
S=zeros(N/2+1,number_of_frame);
for frame=1:number_of_frame,
offset=shift_size*(frame-1);
for n=1:N,
x(n)=s(offset+n)*w(n);
end
X=fft(x,N);
for k=1:N/2+1,
S(k,frame)=20*log10(1+abs(X(k)));
end
end
time=zeros(1,number_of_frame);
for frame=1:number_of_frame,
time(frame)=(window_size/2)/fs+(frame-1)*(shift_size)/fs;
end
frequency=zeros(1,N/2+1);
for k=1:N/2+1,
frequency(k)=(k-1)*fs/N;
end
それと対になるテストスクリプト。
ob$ cat ex.m
clear;
fs=6000;
bits=16;
s=audioread('S.wav');
window_size=512;
shift_size=64;
[S,time,frequency]=spectrogram(s,fs,bits,window_size,shift_size);
imagesc(time,frequency,S);
colormap(flipud(gray));
set(gca,'ytick',[0:500:3000]);
set(gca,'xtick',[0:1:10]);
axis([0 10 0 3000]);
axis xy;
結果が出てくるまで、随分と待たされた。そこで一計。pkg signalを入れちゃえ。やってみたら、哀れな事にgfortranが必要って言われたよ。そんな物は入れたくない。ダメ元で、specgram.mだけを貰ってきた。
octave:12> [D,Fs] = audioread('S.wav');
octave:13> specgram(D,1024,Fs)
これだけで、さっと声紋が出てきた。specgram.mを本体に混ぜてしまえ。
ob$ ls /usr/local/share/octave/4.4.1/m +containers/ gui/ ode/ profiler/ statistics/ @ftp/ help/ optimization/ set/ strings/ audio/ image/ path/ signal/ testfun/ deprecated/ io/ pkg/ sparse/ time/ elfun/ java/ plot/ specfun/ general/ linear-algebra/ polynomial/ special-matrix/ geometry/ miscellaneous/ prefs/ startup/
ここのsignalの下に突っ込んでおいた。単なるスクリプトだから、こういう事を出来るんだろうね。そして、そのソースを見ていたら、声紋に特化した例が出てた。それを少々改造。
ob$ cat mari.m
function mari(sndfile, outfile)
[x, Fs] = audioread(sndfile); # audio file
step = fix(5*Fs/1000); # one spectral slice every 5 ms
window = fix(40*Fs/1000); # 40 ms data window
fftn = 2^nextpow2(window); # next highest power of 2
[S, f, t] = specgram(x, fftn, Fs, window, window-step);
S = abs(S(2:fftn*4000/Fs,:)); # magnitude in range 0<f<=4000 Hz.
S = S/max(S(:)); # normalize magnitude so that max is 0 dB.
S = max(S, 10^(-40/10)); # clip below -40 dB.
S = min(S, 10^(-3/10)); # clip above -3 dB.
colormap('ocean');
H = imagesc (t, f, log(S)); # display in log scale
colorbar();
set (gca, "ydir", "normal"); # put the 'y' direction
if (nargin == 2)
saveas(H, outfile);
endif
endfunction
信号強度をlogにするのはいいんだけど、自然対数とはいかに？ 普通はlog10が馴染のやつだと思うぞ。まあ、barを見なきゃ関係ないけどね。
通常は、音ファイルを読み込んで声紋確認。これはと言うのに出会ったら、記録を残せるようにしたよ。
octave:10> mari('s.wav')
octave:11> mari('s.wav', 'byOctave.png')
byOctave.png 配色は海をイメージしました。
前の方で出てきた、いきなりspecgramで声紋が出て来るやつ、方や、この関数からの返値を受け取る例。この内部の仕組みはどうなってるの？ 現場に当たれ！
specgram.m
function [S_r, f_r, t_r] = specgram(x, n = min(256, length(x)), Fs = 2,
window = hanning(n), overlap = ceil(length(window)/2))
if (nargin < 1 || nargin > 5)
print_usage ();
endif
if (! isnumeric (x) || ! isvector (x))
error ("specgram: X must be a numeric vector");
endif
:
if nargout==0
imagesc(t, f, 20*log10(abs(S)));
set (gca (), "ydir", "normal");
xlabel ("Time")
ylabel ("Frequency")
endif
if nargout>0, S_r = S; endif
if nargout>1, f_r = f; endif
if nargout>2, t_r = t; endif
この関数が呼ばれた時に、C語のmain()宜しく、argc相当だけじゃなくて、（期待する）出力の個数が、nargoutに乗ってやって来るのね。中では、それに合わせて挙動を変えているんだな。
実行スピードが、サンプル実装に比べて随分と速い。その速さの秘密は？
## build matrix of windowed data slices
offset = [ 1 : step : length(x)-win_size ];
S = zeros (n, length(offset));
for i=1:length(offset)
S(1:win_size, i) = x(offset(i):offset(i)+win_size-1) .* window;
endfor
## compute Fourier transform
S = fft (S);
fftすべきデータをwindowで示される窓関数を掛けて事前に用意。それから、一気にfftしてる。forってのは遅いから使うなって典型ですな。
そう言えば、何かする時は行列にデータを用意して、適用ってパターンが多いんで、コマンド列をただ並べたって感じのやつが多いな。
imagesc
上で声紋を書かせると、前回やったpythonと同じ色使いの図が出てきた。どっちかが真似してないかい？ pythonのそれは、見る気が失せるので、octaveのやつを見る。まずはデモデモ。
octave:17> demo imagesc
imagesc example 1:
clf;
colormap ("default");
img = 1 ./ hilb (11);
x = y = -5:5;
subplot (2,2,1);
h = imagesc (x, y, img);
ylabel ("limits = [-5.5, 5.5]");
title ("imagesc (x, y, img)");
:
imagesc example 5: # bug #48879
clf;
img = reshape (1:100, 10, 10);
imagesc (img);
colormap (prism (10));
title ("10 vertical color bars");
どうやらcolormapってので、配色を決定してるようだ。 colormapを調べてみる。
‘colormap (MAP)’ sets the current colormap to MAP. The colormap
should be an N row by 3 column matrix. The columns contain red,
green, and blue intensities respectively. All entries must be
between 0 and 1 inclusive. The new colormap is returned.
Map Description
--------------------------------------------------------------------------
viridis default
jet colormap traversing blue, cyan, green, yellow, red.
：
hot colormap traversing black, red, orange, yellow, white.
cool colormap traversing cyan, purple, magenta.
spring colormap traversing magenta to yellow.
summer colormap traversing green to yellow.
autumn colormap traversing red, orange, yellow.
winter colormap traversing blue to green.
システムで何種類か用意してくれている。春夏秋冬なんてのもあるけど、あちらの人がイメージする季節の色って、日本人とはかけ離れているな。
標準で用意してるマップは64段階になってる。rgbplot で確認できる。
oceanあたりが、夏にふさわしい配色(個人の感覚ですが)と思うぞ。
octave:3> summer(10) ans = 0.00000 0.50000 0.40000 0.11111 0.55556 0.40000 0.22222 0.61111 0.40000 0.33333 0.66667 0.40000 0.44444 0.72222 0.40000 0.55556 0.77778 0.40000 0.66667 0.83333 0.40000 0.77778 0.88889 0.40000 0.88889 0.94444 0.40000 1.00000 1.00000 0.40000
memo
OpenBSDは、sndiodを動かしておくと、音の扱いがまともになる。ただ、ステレオな音源はダメなんだよな。