パワースペクトル計算の2つの方法 (Julia Advent Calendar 2019)

はじめに

これはJulia Advent Calendar 2019の20日目の記事です。「音について何か書きます」と宣言したので、普段Juliaを使って音についてどんなことをしているかを紹介します。今年はパワースペクトルの計算を、(1)信号全体をフーリエ変換する方法と、(2)スペクトログラムを時間平均する方法の2とおりでやってみます。

去年はAdvent Calendar参加記事として以下のようなものを書いていましたので、それからあまり進歩していませんが、ご容赦を。

marui.hatenablog.com

今回のエントリはJulia 1.0.5用に書いています。というのも、Julia 1.2や1.3あたりからLibSndFileなどのインストールがうまくいかない状況になっているからです。（原因はSampledSignals.jlにあるようですが、本稿を執筆している12月19日時点で解決されていません　→2019-12-31 解決しました）

julia> versioninfo()
Julia Version 1.0.5
Commit 3af96bcefc (2019-09-09 19:06 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.6.0)
  CPU: Intel(R) Core(TM) i7-7920HQ CPU @ 3.10GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, skylake)

Julia（とR）で古典的多次元尺度法

昨日、こんなツイートを見た。

以前、県庁所在地と県庁所在地の間の距離が一番近いのどこか気になって、大津と京都が一番近いという無駄知識を得たんだけど、国土地理院さんが「都道府県庁間の距離」というなんのために作成されたのかまったくわからない図表を作成・公開していた。 pic.twitter.com/omv699HW3E

— 森川 真 (@mmww) 2019年11月14日

このデータは以前にも見たことがあって、多次元尺度法のよいサンプルになるかと思っていたのでした。今日はJuliaで古典的多次元尺度法をやってみる。

PortAudio.jlのインストール（2019年11月3日時点）

PortAudio.jlが公式レジストリから消えてしまったのか、これまでのようにインストールしようとすると「見つからない」と言われてしまいます。

(v1.2) pkg> add PortAudio#julia1
ERROR: The following package names could not be resolved:
 * PortAudio (not found in project, manifest or registry)
Please specify by known `name=uuid`.

GitHubにはこれまでどおりコードが置いてありますので、

add https://github.com/JuliaAudio/PortAudio.jl#julia1

とURLを指定することで直接GitHubのレポジトリからインストールすることができます。

AES147 New York技術発表より

10月16日(水)〜10月19日(土)に行われた147th AES Conventionに行ってきました。AES Conventionは大きく製品展示と技術発表に分かれます。僕は技術発表を中心に参加して30件以上の発表を見聞きしましたが、その中で個人的に面白いと思った発表をいくつか書き残しておこうと思います。技術発表の中にもスタジオや放送の現場に近い内容を扱う発表もあれば、より基礎的な研究に焦点を絞った発表もあります。以下、どちらかというと研究方面に偏っています。

carとcdrとcons

学生のときに受けた講義に「プログラミング」がありました。「プログラミングI」「同II」「同III」と半期の科目が3つあり、それぞれ「プログラミングI」は「基本のC言語」、「プログラミングII」は「色々経験するために、関数型言語（LISP）・論理型言語（Prolog）・並列並行処理（AEGISという言語だったような？）」、「プログラミングIII」は「（だいぶ昔のことなのでうろ覚えですが）今後はオブジェクト指向プログラミングが主流になるからC++とJava」というようになっていた気がします。

「プログラミングII」で関数型言語に触れる部分ではLISPを習いました。いま思えばMLとかでも良かったかと思いますが、当時はLISPは実用性のあるメジャーな言語という雰囲気だったのかもしれません。歴史的にも勉強しておいた方がいい言語の一つですしね。

LISPで最初につまづいたのがリスト操作です。いまから思い返すとなぜ理解できないのか分からない、なんて気がしますが、当時はcar、cdr、consでのリスト操作に違和感があり、何が何だか分かりませんでした。おそらくcar+cdrとconsが対称な操作であることが分かっていなかったのではないかと思います。

たとえば

(car '(1 2 3))   -> 1
(cdr '(1 2 3))   -> (2 3)
(cons 1 '(2 3))   -> (1 2 3)

について、carはリストの最初の要素をとってくる、cdrはリストの残りをリストとして得る、consは引数2の先頭に引数1を挿入するという動作をするので、当時の自分の理解を図解すると以下のような感じです。

この簡単なことが当時の僕にはよく分からず、

(cons '(1 2) '(3 4))   -> (1 2 3 4)

を期待しているのに

(cons '(1 2) '(3 4))   -> ((1 2) 3 4)

となるので「先頭に値を入れてくれるのではないのか？」などと混乱していました。同じ頃にMatlabも使い始めていたので、

>> A = [1 2];
>> B = [3 4];
>> [A B]
ans =
     1     2     3     4

と同じような動作を期待していたのだと思います。

LISPに触れるのは半期の授業の1/3だけなので4回程度。話題はあっという間に他の言語に移ってしまったので「論理型言語おもしろいけど、家系図つくって実用性あるのか？」だの「並列・並行は意味が分からない、円卓の哲学者は馬鹿か？」などと言い始め、期末試験が終わるとそれら全てを忘れてしまい、卒業までずっとC言語を使っていました。

その後、何年も経ってからLISPに再入門すると、car+cdrとconsの見え方が変わっていました。値やリストを抽出・挿入という考え方から、括弧の位置が動くという考え方になりました。図解すると以下のような感じです。

carやcdrは括弧が1つ右に動く感覚です。たとえば(car '(1 2 3))は(1 2 3)の先頭の括弧が1 (2 3)と1つぶん右に移動し、そのときに括弧外に出たものがcarで、括弧内に残っているものがcdrです。一方でconsは括弧が1つ左に動く感覚です。リストは1つのまとまりなので、(cons '(1 2) '(3 4))をすると、(3 4)の先頭の括弧が(1 2)のまとまりを超えて左側に動き((1 2) 3 4)となるのでした。

結局、(1 2)と(3 4)から僕が期待したとおりの(1 2 3 4)を作るためには、(1 2)の各要素を取り出しつつconsを使い

(cons (car '(1 2)) (cons (car (cdr '(1 2))) '(3 4)))   -> (1 2 3 4)

とすることになります。（これをするために(append '(1 2) '(3 4))が用意されています）

アダマール行列（N = 2^kのときのみ）

StautnerとPucketteが提案したFeedback Delay Networkという人工リバーブがあります。このFDNを実装するときに必要になるものに、N×Nの行列があります。これを行列Gとしたとき、各要素g_{i,j}はディレイラインiからディレイラインjへのフィードバック・ゲインが入っています。この行列はいろいろと選択肢があるのですが、J. Smithによると「アダマール行列も使えるよ」ということでした。

MaxやPureDataでFDNを実装するときにアダマール行列が欲しかったので、Juliaで計算したものを入力して使おうと思いました。ただ、MatlabやOctaveにあるようなhadamard関数がJuliaにはなかったので、作りました。参考にしたのはWikipediaに載っていた“他の生成法”。MatlabだとN=2^k * pでp=1, 12, 20, 28のときの計算ができるけど、とりあえず今回はp=1のときだけ。

"""
    hadamard(n)

Compute the Hadamard matrix of size `n`-by-`n`.
`n` must be a positive integer of power of 2.
"""
function hadamard(n)
  function hadamardF(n)
    if n <= 1
      return [0 1];
    else
      return [zeros(1, 2^(n-1))   ones(1, 2^(n-1))
              hadamardF(n-1)      hadamardF(n-1)];
    end
  end
  
  F = hadamardF(floor(Int, log2(n)));
  FF = F' * F;
  FF[(FF .% 2) .== 1] .= -1;
  FF[(FF .% 2) .== 0] .= 1;
  return map(Int, FF);
end

以下のように使えます。

julia> hadamard(4)
4×4 Array{Int64,2}:
 1   1   1   1
 1  -1   1  -1
 1   1  -1  -1
 1  -1  -1   1

ただ、アダマール行列を計算するプログラムはすでに色々とあります。上記のコードよりもっと高機能なものが以下のリンク先にあるので、それを使うのもよいですね。

nemocas.org

github.com

線形予測符号（LPC）でスペクトル包絡を計算する

Juliaを使った音の処理については、これまでにいくつかのエントリを書きましたが、今日は線形予測符号（LPC; Linear Predictive Coding）を用いてスペクトル包絡を求めてみます。ここでは線形予測符号の説明は省きますが、音声信号処理の分野では昔からよく使われている方法で、Juliaコードを書くにあたっては以下のページを参考にしています。

• 『Miyazawa's Pukiwiki』より「MATLAB Note/音声の分析」（http://speechresearch.fiw-web.net/73.html）
• 奈良先端大学の講義「音情報処理論」第2回「音声の特徴抽出」資料（https://ahclab.naist.jp/lecture/2017/sp/material/slide.pdf）
• 『人工知能に関する断創録』より「線形予測分析（LPC）」（http://aidiary.hatenablog.com/entry/20120415/1334458954）

なお、JuliaでLPCを自作することはせず、DSPパッケージに入っているlpc()をありがたく利用します。

準備

まず、使用するパッケージ群の読み込みを行います。

using FileIO: load
import LibSndFile
using DSP
using FFTW
using Plots

音ファイルの読み込み

LibSndFileを使って音ファイルを読み込んできます。今回は音声ではなく、とある管楽器の音を使ってみます。

snd = load(expanduser("~/instrument.wav"))

x = snd.data[1:60000];
fs = snd.samplerate;
t = [0:length(x)-1] / fs;

plot(t, x, ylim=[-1, 1],
     xlabel="Time (s)", ylabel="Amplitude", legend=false)

楽器音を読み込んで、そこから音信号xと標本化周波数fsを取り出します。また波形のプロットに使いたいので、xの長さに合わせた各サンプルの時刻tを計算します。plot()で波形表示をしています。

音声の場合はフォルマントを把握しやすくするために高域を強調する「プリエンファシス・フィルタ」をかけるそうなのですが、今回は楽器音を使っているので、そのままにしておきます。

パワースペクトルの計算

パワースペクトルの計算は必須ではありませんが、スペクトル形状を確認したいので計算しています。

y = fft(x) / length(x);
yy = abs.(y);
yy = yy[1 : div(length(y), 2)+1];
yy[2:length(yy)-1] .*= 2;
f = range(0, fs/2, length=length(yy));

plot(f, 20*log10.(yy),
     alpha=0.5, label="Power Spectrum",
     xlabel="Frequency (Hz)", ylabel="Power (dB)")

フーリエ変換（fft）をして、正の周波数領域の振幅スペクトルだけを抜き出す処理をしています。

各行をなぞっていきます。まずabs()で振幅スペクトルを得て（2行目）、DCからナイキストまでの成分を抜き出します（3行目）。すると負の周波数側の成分を捨てたことになるので2倍して補正します（4行目）。yyは各周波数ビンの振幅が入っているので、各ビンの中心周波数を求めます（5行目）。

プロットでは、あとでLPCの結果を重ねたいので、alphaを指定して振幅スペクトルの色を薄くしています。

LPC

ここからが本番です。DSPパッケージにlpc()が入っているので、それを使います。LPCの次数は適当に決めて、20次でやってみます。lpc()を使うとフィルタ係数aが求められます。

p = 20
a, e = lpc(float(x), p);
env_x = filt([1], [1; a], [1; zeros(length(x)-1)]);
env_y = fft(env_x) / length(env_x);
env_yy = abs.(env_y);
env_yy = env_yy[1 : div(length(env_y), 2)+1];
env_yy .*= 2;

plot!(f, 20*log10.(env_yy),
      width=1.5, label="Spectral Envelope (LPC p=$(p))")

フィルタ係数からエンベロープの形状を計算するのにDSPパッケージのfreqz()を使えることをすっかり忘れていたので、上記ではフィルタにインパルスをつっこんで周波数応答を求めています。

おわりに

他にもケプストラムを使う方法もありますが、そちらについてはまた別の機会に試してみたいと思います。

上記の方法は厳密に正しくない部分があるかもしれません（たとえばFFT後に正負のエネルギを合算するところは本当に必要なのか分かりません）。「ここが間違ってるよ」「こうした方がいいよ」といった意見ありましたら、コメント欄からでも連絡ください。