線形予測符号（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()で波形表示をしています。

f:id:amarui:20190814135831p:plain

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

パワースペクトルの計算

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

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を指定して振幅スペクトルの色を薄くしています。

f:id:amarui:20190814140424p:plain

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))")