MatlabのfreqzもどきをJuliaで

MatlabOctaveにはfreqzという関数があり、FIR/IIRフィルタの周波数特性・位相特性を計算しグラフ表示してくれます。PythonでもSciPy.signalに同名の関数があります。似たことをJuliaでやってみました。

今回は双二次フィルタを見てみましょう。標本化周波数44,100 Hz、中心周波数1,000 Hz、ゲイン12 dB、Q=2.0のピークフィルタを準備しました。

samprate = 44100;
B = [1.1417, -1.9797, 0.8583];
A = [1.0356, -1.9797, 0.9644];

特性を見たい周波数範囲を設定し、それをzさらにsへと変換して使います。ここでは16 Hz〜ナイキスト周波数を256段階の対数分割してみました。

# calculate frequency
# (using logspace since we are going to use semilogx for plotting)
fL = 16;
fH = samprate/2;
f = logspace(log10(fL), log10(fH), 256);
z = f / sampfreq * 2*pi;
s = exp.(-1.0im * z);   # convert from digital frequency

次に、各周波数での伝達関数の値を計算します。下の式に代入して計算するだけ。

f:id:amarui:20170927203946p:plain
[tex: H=\frac{b_0 + b_1 s + b_2 s2}{a_0 + a_1 s + a_2 s2}]

# calculate transfer function
H = (B[1] + B[2].*s + B[3].*s.*s) ./ (A[1] + A[2].*s + A[3].*s.*s);

Polynomialパッケージを使って汎用性を高めても良いのですが、僕の考えている使用用途は双二次フィルタだけなので、このままで。

さて、振幅特性をグラフにしてみましょう。abs.()で絶対値を計算してからデシベルにします。

using PyPlot
semilogx(f, 20*log10.(abs.(H)));
grid();
xlabel("Frequency (Hz)");
ylabel("Gain (dB)");
savefig("freqz_abs.png")

f:id:amarui:20170927145542p:plain

次に位相特性もグラフにしてみましょう。今度はangle.()で位相を抜き出し、ラジアンを角度に変換しています。

semilogx(f, angle.(H)*180/pi);
grid();
xlabel("Frequency (Hz)");
ylabel("Phase (angle)");
savefig("freqz_phase.png")

f:id:amarui:20170928163008p:plain

Matlabfreqzのグラフと見比べてみましたが、だいたい同じ結果になっているようでした。Matlab版は横軸が線形になっているので見た目が違いますが。 f:id:amarui:20170928162206p:plain

振幅スペクトルのリアルタイム可視化

とある場所で振幅スペクトルのリアルタイム可視化をデモすることになったので、久しぶりにProcessingのプログラムを書きました。Minimに入っていた例をほぼそのまま使っています。(左右方向は20 Hz〜ナイキストを対数軸で、上下方向は-60〜0 dBを表示しています)

f:id:amarui:20170919152525p:plain

import ddf.minim.analysis.*;
import ddf.minim.*;

Minim minim;  
AudioPlayer jingle;
FFT fft;

void setup() {
  size(1024, 256, P3D);
  
  minim = new Minim(this);
  jingle = minim.loadFile("hogehoge.mp3", 1024);
  jingle.loop();
  
  fft = new FFT(jingle.bufferSize(), jingle.sampleRate());
}

void draw() {
  // set background color
  background(#0D2B36);  

  // compute FFT
  fft.forward(jingle.mix);
  
  // draw foreground
  stroke(#EEE8D5);
  strokeWeight(2);
  for(int i = 0; i < fft.specSize(); i++) {
    float y = fft.getBand(i) * 2.0 / fft.specSize();
    float xlog = log(float(i) / fft.specSize()) / log(10);
    float ylog = 20 * log(y) / log(10);

    // adjust x/y positions to fit into window
    xlog = (3.0 + xlog) / 3.0 * width;
    ylog = height - (ylog+60)/60 * height;

    // color based on amplitude
    stroke(color(255*ylog/height, 255*ylog/height*0.8, (255-255*ylog/height)*0.5));

    line(xlog, height, xlog, ylog);
  }
}

どれも1行〜2行の変更だけでできるので、楽しくてついつい必要以上のことをやってしまう。アレンジしだいで色々できるし、子どもの頃にこういうのがあったら無限にやってただろうな。(下のスクリーンショットは、周波数と振幅を極座標系で表示し、左右対称にコピーしたもの。蝶みたいに見える。)

f:id:amarui:20170919155040p:plain

本日の給油

給油日 オドメーター (km) 給油量 (L) 単価 (円/L) 燃費 (km/L) 距離単価 (円/km)
2017-03-18 8559.6 2.76 125.72 45.94 2.74
2017-04-29 8703.9 2.77 127.80 52.09 2.45
2017-04-29 8905.9 3.21 128.97 62.93 2.05
2017-04-30 9032.5 2.20 131.82 57.55 2.29
2017-05-01 9190.3 2.70 127.04 58.44 2.17
2017-05-01 9325.7 2.52 128.97 53.73 2.40
2017-05-01 9469.7 2.17 129.03 66.36 1.94
2017-05-02 9625.3 3.71 133.96 41.94 3.19
2017-05-02 9787.8 2.81 125.98 57.83 2.18
2017-09-15 10004.7 3.09 123.96 53.03 2.34

10,000 km乗りました。「3ヶ月で!!」とかならかっこいいけど、2,000 km/年なので自慢できるものではないですが……。

MacBook Pro (15in, 2017)の3 chマイク

MacBook Pro (15in, 2017)ですが、「Audio MIDI設定」の画面でマイク設定を確認してみると4 ch入力ができるかのように書いてあります。

f:id:amarui:20170903093824p:plain

最高で96 kHz / 24 bit / 4 chまでいけるとなれば、マイクロホンとヘッドアンプの性能によっては本体だけでもかなり面白いことができるのでは、と思います。(実際、電気グルーヴGarage Bandと内蔵マイクで多くの部分を録音したアルバムなんてのもありますし)

カタログスペックを調べてみると、

  • Stereo speakers with high dynamic range(広ダイナミックレンジのステレオスピーカ)
  • Three microphones(マイクロホン×3)
  • 3.5 mm headphone jack(3.5 mmヘッドホンジャック)

となっていました。「マイクロホン×3」って! 奇数のADCを載せるのは面倒(かえって高くつく?)から、マイクは3つだけどDACは4 chのものを搭載しておこうということなのでしょうか。

試しにAudacityを使って4 ch録音をしてみました。画像では見にくいかもしれませんが、たしかに1〜3 chは録音された波形が違い、4 ch目は3 ch目と同じ波形です。

f:id:amarui:20170903095622p:plain

これら3つのマイクを何に使っているかですが、「システム環境設定」→「サウンド」でも表示されるように(「Audio MIDI設定」でマイクを2 chにしたときのみ表示)、これらを使って暗騒音や環境音の除去をやっているんでしょう。話者方向の推定とかにも使えるかもですね。

MacBookから映像が出せなくて困ったはなし

Apple MacBook Pro (2017)とエレコムのDST-C01SVとで講義デビュー。

DST-C01SVは、ギガビット・イーサ、VGAHDMI、Mini Display Port、ステレオミニジャック(入出力対応)、USB 3.0×3口、SDスロット、microSDスロットに加えて、充電専用のUSB Type-C端子までついているという、現時点で敵なしのすごいやつです。MacBook AirのようにType-Cポートが少ない機種でも充電しながら他のものも接続できて、これ一つでどこでも行ける。欲を言えば、本体から出てるType-Cケーブルは(Apple純正の電源アダプタのように断線防止のために)着脱できるようにしておいて欲しかった。

www2.elecom.co.jp

ただ今日は、映像を使おうとすると、HDMIVGAも認識しようとして固まった後、15秒ほどして認識を諦めてしまうという現象に遭遇。いろいろ試すも対処法が見つからず結局ダメ。

そのためなんと今日の講義は3時間以上にわたって「実物投影機でMacBookの液晶画面を撮影しながらそれをプロジェクション」というすごいことになってしまった。RetinaディスプレイフルHDのプロジェクタだったので画質はまぁまぁでしたが、講義室の機材卓では実物投影機からの映像入力とPCからのオーディオ入力が同時利用できなかったので、絵と音を一緒に見せられない。音の波形を拡大して一部分だけループ再生、みたいなことができませんでした。

帰り着いて、なんとかならないかとググったら以下のページが。

digmusik.blogspot.com

……結局、僕の場合は「NVRAMリセット(PRAMリセット)」でHDMIに映像が出力できました。明日からの3日間は(その翌日からの別の3日間出張も)なんとかなりそうです。

ただし、いったん接続状態で本体がスリープしてしまうと再接続ができず、またNVRAMリセットが必要になる気がする。完全解決はまだですなー。おそらくmacOSのアップデートが必要です。

(後日談:OSアップデートで解決しました。今は快適です。)

同名の中華香辛料が二種類あった?

僕が好きな辛味調味料のうち、香りの高さで選んでいるものが老干媽(ラオガンマ)が作っている香辣脆油辣椒(シャンラーツィユーラージャオ)です。定期的に横浜中華街まで買いに行っていることは下記のエントリにも書きました。

marui.hatenablog.com

東京の東側に住んでいるので横浜中華街まで行こうとすると電車賃も時間もかかります。近所で買えればいいのにと思っていたところ、3年ほど前から業務スーパーで油辣椒と風味豆豉油制辣椒が扱われるようになりました。香辣脆油辣椒ほどではないものの、けっこう好きな調味料なので定期的に購入していました。

そして、半年前くらいからは業務スーパーで香辣脆油辣椒も買えるようになりました!! 老干媽の製品の中でも唐辛子の辛味(辣)より山椒の痺れ(麻)が効いていて、大好物の調味料です。瓶のシールはデザインが変わったようですが、商品名は変わらず、おなじみの顔写真もついています。

……で、自宅で使ってみると、なんだか味も香りも違います。辣も麻も弱く、油っぽさとタマネギの香りが強く感じられます。香ばしいタマネギっぽさが向上しているので、好きな人にはたまらないのでしょうが、個人的にはもっとガッツリとした辣と麻の香りの強さが好みですし、以前と味が変わってしまったのを残念に思いました。

そして先日、ちょうど買い置きの香辣脆油辣椒がなくなったこともあり、久しぶりに横浜中華街に行くことにしました。今回も栄興號でひとつあたり270円の香辣脆油辣椒を12本購入。だいたい350円で売っている他店よりも安いのにていねいな対応と梱包で、大好きな店です。

f:id:amarui:20170818171744j:plain

帰宅してさっそく食べ比べてみました。やはり味が違いますし、ラベルも違います。業務スーパーで買ってきたものは神戸物産が輸入者で、中華街で買ってきたものは友盛貿易が輸入しているものです。この1瓶だけでなく6瓶ほどが同じ傾向だったので、製造者が輸入者のリクエストに応じて原料の配合を変えているのかもしれません。インスタントのカップうどんなどでも東日本と西日本で出汁の味を変えている例がありますし、ラベルを見たら区別できるので、自分の好みに合わせて選べる状況は悪いことではないのかと(僕は友盛貿易扱いのシャープな味が好みなので中華街に通います)。

音ファイルの読み込みと単純な可視化

ここ数年の僕はプログラミングというと音と統計ばかりやっているのですが、特に音の分析や加工に関してよく使うものはMatlabPython、Rで実装したので、次はCommon LispかJuliaで実装してみようかと思います。最終的にはCommon Lispでなんでもやるのが夢なのですが、今回はJuliaで。

基本的なことがら

音ファイルを扱えないと、その後の分析も何もできませんので、いの一番にやることはJuliaでWAVファイルの読み込みです。音に関するパッケージをPkg Serverで探してみると、使えそうなパッケージはいくつかあります。その中でもLibSndFileを使えるようにしてくれるパッケージが(現時点で丸井綜研が使っている)Julia 0.6のテストを通過してるので良さそうです。

インストール

インストールは簡単で、Juliaのコマンドラインから

Pkg.add("LibSndFile")

とやるだけ。初回に一度やれば、次からはこの部分は不要です。

使用準備

using LibSndFile

でパッケージの読み込みができます。初回はプリコンパイルなどで少し時間がかかるかもしれません。

音の読み込み

読み込みにはloadを使います。JuliaのコマンドラインREPLではいろいろな情報を表示してくれます。以下ではCDから抜粋した15秒ほどのWAVファイルを読み込んでいますが、1チャンネルあたり587,264サンプルあり、それが2チャンネルあること、全体で約13.3秒、サンプリング周波数は44,100Hzなどの情報に加えて、波形の絶対値?をプリティプリントまでしてくれます。

julia> x = load("gershwin_piano_concerto.wav")
587264-frame, 2-channel SampleBuf{FixedPointNumbers.Fixed{Int16,15}, 2}
13.316643990929705s sampled at 44100.0Hz
▅▆▆▆▆▇▆▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▆▆▇▇▇▇▇▇▇▇▇▇▇▇▇▆▆▇▆▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▆▇▇▇▇▇▇▇▇▇▇▇▆▇▆▆▇▆▆▆▆▅▄
▅▆▆▆▆▆▆▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▆▇▇▆▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▆▆▆▆▆▆▆▅▄

IJuliaノートブック上では、なんとグラフィカルな波形表示と再生機能までしてくれます。(下はそのスクリーンショットです)

f:id:amarui:20170810150021p:plain

格納先の変数はSampledSignals.SampleBuf型で、datasamplerateの二つのフィールドを持っています。

julia> x.data   # 各サンプル値(2chファイルなので2列です)
587264×2 Array{FixedPointNumbers.Fixed{Int16,15},2}:                            
 -3.0e-5Q0f15    0.0Q0f15                                                       
  0.0Q0f15       3.0e-5Q0f15                                                    
 -3.0e-5Q0f15    3.0e-5Q0f15                                                    
 -3.0e-5Q0f15    3.0e-5Q0f15                                                    
 -3.0e-5Q0f15    6.0e-5Q0f15                                                    
  ⋮                 
  6.0e-5Q0f15   -6.0e-5Q0f15            
  0.0Q0f15       3.0e-5Q0f15            
 -6.0e-5Q0f15    3.0e-5Q0f15            
 -6.0e-5Q0f15    6.0e-5Q0f15            
 -6.0e-5Q0f15    6.0e-5Q0f15            
 -6.0e-5Q0f15    3.0e-5Q0f15            

julia> x.samperate   # サンプリング周波数
44100.0

波形の表示

音の加工などをする前に波形を見てみましょう。

x1 = x.data[:, 1];   # チャンネル1のみ抜き出し
t = (0 : length(x1)-1) / x.samplerate;   # 各サンプルの時刻を求める

using PyPlot
plot(t, x1);
ylim(-1, +1);
grid();
title("Music Excerpt")
xlabel("Time (s)")
ylabel("Amplitude")

f:id:amarui:20170810164125p:plain

ひとまずはここまで。Octaveを使っているのとそれほど変わりない感覚です。

応用例(周波数分析)

おまけとして、リコーダーの音の倍音構造を見てみようと思います。変数名がぐちゃぐちゃだけど許して。

y = load("recorder.flac");   # FLACも読み込める
fs = y.samplerate;
y1 = y[1:Int(fs*1.8)];   # 冒頭から1.8秒分だけ使う(約524Hzの吹奏音が入っています)

y2 = zeros(nextpow(2, length(y1)));   # y1より長い、2の累乗の長さの零ベクトルを作る
y2[1:length(y1)] = y1;   # その先頭にy1を入れる(つまりzero-padding)
Y2 = abs.(fft(y2));   # フーリエ変換して絶対値を計算
Y2 = (Y2 / length(y2)) .^ 2;
Y3 = Y2[1 : Int(length(Y2)/2) + 1];
Y3[2:end-1] = Y3[2:end-1] * 2.0;   # 正の周波数成分と負の周波数成分をまとめる
Y4 = 10*log10.(Y3);   # デシベルに変換

f = linspace(0, y.samplerate/2, length(Y3));   # 各ビンの周波数を計算

plot(f, Y4-maximum(Y4));   # 最大値が0 dBになるようにしてグラフ描画
xlim(0, 15000);   # X軸・Y軸の範囲指定
ylim(-80, +3);
grid();   # グリッドを描画
title("Power Spectrum");   # タイトルや軸のラベルを付ける
xlabel("Frequency (Hz)");
ylabel("Power (dB)");

savefig("output3.png");   # PNGファイルへの保存

f:id:amarui:20170810173526p:plain

この音は第2・第4倍音が少なめですね。