前回実装した、NNLSのLawson–Hansonアルゴリズムで、能動集合（passive set）に制限した最小二乗問題を解くときにコレスキー分解を行っている。
辞書行列Eが、各半音ごとに独立した倍音スペクトル分布を列として持つため、ほぼフルランク（すべて線形独立）であることを前提にしている。
コレスキー分解とは、対称正定値行列 を次のように分解する方法である。
ここで：
つまり、複雑な行列を「三角形の形」に分解できる、というのがコレスキー分解のポイントである。
行列の連立方程式を解くとき、そのままでは計算が重かったり不安定だったりする。
しかし三角行列は「順番に代入していくだけ」で解ける。
例として、次の行列を考える：
この行列は対称かつ正定値である。
コレスキー分解すると、
となり、
が確認できる。
前回Pythonで実装したNNLSをWindowsプログラムに組み込んで、リアルタイムで処理できるか検証する。
Pythonでは、scipy.optimize.nnlsを使用したが、C++に組み込むため、Lawson–Hansonのアルゴリズムを独自実装した(GPT-5で実装)。
static void nnls_lawson_hanson_colmajor( const double* A, int m, int n, // A: m×n（列優先）A[j*m + i] = A_ij const double* b, // b: m std::vector<double>& x, // x: n（出力） std::vector<double>& workAx, // Ax: m std::vector<double>& workRes, // r = b - A x: m std::vector<double>& workW, // w = A^T r: n std::vector<int>& passive, // 1=Passive(P), 0=Active(Z) std::vector<int>& active // 1=Active(Z), 0=Passive(P) ) { const int MAX_IT = 3 * n + 5; const double TOL = 1e-10; x.assign(n, 0.0); passive.assign(n, 0); active.assign(n, 1); // w = A^T b workAx.assign(m, 0.0); workRes.assign(m, 0.0); workW.assign(n, 0.0); for (int j = 0; j < n; ++j) { double dot = 0.0; const double* col = A + j * m; for (int i = 0; i < m; ++i) dot += col[i] * b[i]; workW[j] = dot; } int it = 0; while (it++ < MAX_IT) { // Z（active）から最大正の w を持つ列を選択 double wmax = 0.0; int t = -1; for (int j = 0; j < n; ++j) if (active[j] && workW[j] > wmax) { wmax = workW[j]; t = j; } if (t < 0 || wmax <= TOL) break; // 収束 passive[t] = 1; active[t] = 0; // 反復内ループ：負成分が出たら P→Z に戻す std::vector<int> Pidx; Pidx.reserve(n); for (int j = 0; j < n; ++j) if (passive[j]) Pidx.push_back(j); // 最小二乗（正規方程式）を P のみで解く std::vector<double> G, g; auto solve_on_P = [&](std::vector<double>& xP) -> bool { const int p = (int)Pidx.size(); if (p == 0) return false; G.assign(p * p, 0.0); g.assign(p, 0.0); // G = A_P^T A_P, g = A_P^T b for (int a = 0; a < p; ++a) { int ja = Pidx[a]; const double* colA = A + ja * m; for (int b2 = a; b2 < p; ++b2) { int jb = Pidx[b2]; const double* colB = A + jb * m; double s = 0.0; for (int i = 0; i < m; ++i) s += colA[i] * colB[i]; G[a * p + b2] = G[b2 * p + a] = s; } double sb = 0.0; for (int i = 0; i < m; ++i) sb += colA[i] * b[i]; g[a] = sb; } // xP ← 解 xP = g; if (!cholesky_solve_spd(G, p, xP)) return false; return true; }; // 旧解（P成分） std::vector<double> x_old = x; for (;;) { // 解候補 std::vector<double> xP; if (!solve_on_P(xP)) break; // 候補を全成分ベクトルへ詰める std::vector<double> xCand(n, 0.0); for (int k = 0; k < (int)Pidx.size(); ++k) xCand[Pidx[k]] = xP[k]; // すべて非負なら採用 bool allPos = true; for (int k = 0; k < (int)Pidx.size(); ++k) if (xCand[Pidx[k]] <= 0.0) { allPos = false; break; } if (allPos) { x.swap(xCand); break; } // α を求めて0に落ちる変数を Z へ戻す double alpha = 1.0; for (int k = 0; k < (int)Pidx.size(); ++k) { int j = Pidx[k]; if (xCand[j] <= 0.0) { double denom = xCand[j] - x_old[j]; if (denom != 0.0) { double a = x_old[j] / (x_old[j] - xCand[j]); if (a > 0.0 && a < alpha) alpha = a; } } } for (int j = 0; j < n; ++j) x[j] = x_old[j] + alpha * (xCand[j] - x_old[j]); // 0 になったものを Z へ for (int k = 0; k < (int)Pidx.size(); ++k) { int j = Pidx[k]; if (x[j] <= 1e-16) { x[j] = 0.0; passive[j] = 0; active[j] = 1; } } // P の再構成 x_old = x; Pidx.clear(); for (int j = 0; j < n; ++j) if (passive[j]) Pidx.push_back(j); if (Pidx.empty()) break; } // w = A^T (b - A x) std::fill(workAx.begin(), workAx.end(), 0.0); for (int j = 0; j < n; ++j) if (x[j] != 0.0) { const double* col = A + j * m; double cj = x[j]; for (int i = 0; i < m; ++i) workAx[i] += col[i] * cj; } for (int i = 0; i < m; ++i) workRes[i] = b[i] - workAx[i]; for (int j = 0; j < n; ++j) { if (!active[j]) { workW[j] = 0.0; continue; } const double* col = A + j * m; double dot = 0.0; for (int i = 0; i < m; ++i) dot += col[i] * workRes[i]; workW[j] = dot; } } }
リアルタイムでも問題なく処理できることがわかった。
Debugビルドしたバイナリでも遅延なく処理できる。
スペクトログラム解析VST3プラグインに和音の解析機能を加えたい。
ミックスされた音源から和音の構成音を取り出すのは容易ではなく、古典的にはCQT や HPCPをピッチクラス化して推定する方法がある。
この方法では、倍音が構成音と誤検知される場合があり、精度を上げるのが難しい。
最新の精度の高い手法は、Deep Chromaのようなディープラーニングベースの手法が用いられる。
ディープラーニングベースの手法は計算量が大きくスマホでのリアルタイム解析には向かないため、軽量である程度の精度がでる手法がないか調査したところ、「Approximate Note Transcription for the Improved Identification of Difficult Chords」という方法を見つけた。
上記の論文の手法では、NNLS Chromaという特徴量からコードを推定する。
NNLSは、Non-Negative Least Squares（非負最小二乗法）の略で、解が正の値のみという制約付き最適化問題である。
CQTのスペクトルから、A0からG#6までのどの音が鳴っているかを最適化問題として求める。
スペクトルは複数の音が混ざったものとして、理想の倍音構成の各音が、それぞれどの強さで鳴っているかを求める。
NNLSの解は、解析的には求めることはできないため、反復法で解く必要がある。
NNLSを解くアルゴリズムには、Solving Least Squares Problems(Lawson, C.L. & Hanson, R.J.)や、その改良アルゴリズムA Fast Non‐Negativity‐Constrained Least Squares Algorithm(Bro, R. & de Jong, S.)が用いられる。
後者は、scipy.optimize.nnlsに実装されている。
CQTからピッチクラス化した場合と、NNLS Chromaを使用した場合の比較を行う。
まずは、CQTからピッチクラス化した結果を示す。
音源は、librosa付属のnutcracker(くるみ割り人形「金平糖の精の踊り」)を使用する。
この曲のキーはEマイナーである。
import librosa file_path = librosa.ex("nutcracker") target_sr = 11025 y, sr = librosa.load(file_path, sr=target_sr, mono=True) C = librosa.cqt( y=y, sr=sr, hop_length=2048, bins_per_octave=36, n_bins=36 * 7, ) librosa_chroma = librosa.feature.chroma_cqt( C=C, sr=sr, hop_length=2048 ) librosa.display.specshow( librosa_chroma, sr=sr, hop_length=2048, x_axis="time", y_axis="chroma", cmap="viridis", )
結果
次に、NNLS Chromaからピッチクラス化した結果を示す。
NNLSの実装は、論文にならいChromaとBass chromaをそれぞれ求める。
@dataclass class NNLSChroma: sr: int = 11025 hop_length: int = 2048 bpo: int = 36 n_octaves: int = 7 s_param: float = 0.6 preprocessing: str = "std" chroma_norm: int = 3 def compute(self, y): n_bins = self.bpo * self.n_octaves fmin = librosa.midi_to_hz(21) tuning = librosa.estimate_tuning(y=y, sr=self.sr, bins_per_octave=self.bpo) C = librosa.cqt( y=y, sr=self.sr, hop_length=self.hop_length, fmin=fmin, n_bins=n_bins, bins_per_octave=self.bpo, tuning=tuning, ) Y = np.abs(C).T half_span = 18 proc_Y = np.zeros_like(Y) for t in range(Y.shape[0]): mu, std = _running_mean_std_1d(Y[t], half_span) sub = Y[t] - mu sub[sub < 0] = 0.0 if self.preprocessing == "std": std_safe = np.where(std > 1e-6, std, 1.0) proc_Y[t] = sub / std_safe else: # "o" or "sub" proc_Y[t] = Y[t] if self.preprocessing == "o" else sub E = _build_note_dictionary( n_bins, bpo=self.bpo, n_octaves=self.n_octaves, s=self.s_param, fmin=fmin ) n_semitones = self.n_octaves * 12 S = np.zeros((proc_Y.shape[0], n_semitones), dtype=float) for t in range(proc_Y.shape[0]): x, _ = nnls(E, proc_Y[t]) S[t] = x chroma = np.zeros((S.shape[0], 12), dtype=float) bass_chroma = np.zeros((S.shape[0], 12), dtype=float) pc_offset = 21 % 12 # fmin が A0 のため 9 for t in range(S.shape[0]): for i in range(n_semitones): chroma[t, (i + pc_offset) % 12] += S[t, i] * TREBLE_WINDOW[i] bass_chroma[t, (i + pc_offset) % 12] += S[t, i] * BASS_WINDOW[i] for v in (chroma[t], bass_chroma[t]): if self.chroma_norm == 1: m = v.max() if m > 1e-6: v /= m elif self.chroma_norm == 2: s_norm = v.sum() if s_norm > 1e-6: v /= s_norm elif self.chroma_norm == 3: n = np.linalg.norm(v) if n > 1e-6: v /= n times = librosa.frames_to_time( np.arange(S.shape[0]), sr=self.sr, hop_length=self.hop_length ) return times, chroma, bass_chroma, S, tuning
extractor = NNLSChroma(
sr=target_sr, preprocessing="std", s_param=0.6, chroma_norm=3
)
times, nnls_chroma, nnls_bass_chroma, semitone, tuning = extractor.compute(y)
librosa.display.specshow(
nnls_chroma.T,
sr=sr,
hop_length=extractor.hop_length,
x_axis="time",
y_axis="chroma",
ax=axes[1],
cmap="viridis",
)
librosa.display.specshow(
nnls_bass_chroma.T,
sr=sr,
hop_length=extractor.hop_length,
x_axis="time",
y_axis="chroma",
ax=axes[2],
cmap="viridis",
)
結果
※比較のため、CQTの結果も並べている。
NNLS Chromaの結果は、CQTと比べてピッチクラスがすっきりしており、構成音が分かりやすくなっている。
また、ピッチクラス化前のスペクトル( (d)と(e) )を比べると、(e)のCQTの結果では高域が倍音により密になっているのに比べて、(d)のNNLS Chromaでは低域から高域まで疎な音で構成されているのが分かる。
ミックスされた曲の構成音を推定する手法について調査した。
比較的軽量で精度が高そうな手法であるNNLS Chromaを実装してみた。
CQTからそのままピッチクラス化するよりも構成音が分かりやすくなることが確認できた。
NNLSの解を反復法で求めるにはO(m * n^2)※(m=ビン数252, n=音符数84)の計算量がかかるため、リアルタイムで処理できるか検証したい。
先日、実装したリアルタイムスペクトログラム解析VST3プラグインをPayhipで公開した。
Audio Spectrogram Monitor は、トラックに挿入するだけで使えるシンプルかつゼロレイテンシーのスペクトログラムモニターです。 72分割/オクターブの定Q変換(CQT)ディスプレイに、鍵盤の白黒を模したピアノロール軸、DAWと同期する小節/拍ガイドライン、オプションのピッチライン、さらにMIDIノートオーバーレイを備えています。耳コピ（採譜）や演奏練習、サウンドデザイン、ミックス分析に最適です。
直感的なUI操作
Shift + ホイール: 時間方向のズームCtrl + ホイール: ピッチ方向のズームホイール: 周波数範囲のスクロールCtrl + 0: ズーム/スクロールをリセットC:\Program Files\Common Files\VST3自作スペクトログラム解析VST3プラグインにマウスホイールによる拡大縮小とスクロール処理を実装した。
はじめ解析範囲をVSTのパラメータに使用としたが、マウスホイールの方が操作しやすいのでそちらで実装することにした。
解析範囲を動的に変更するのは処理が複雑になるため、表示範囲のみ変えることで対応した。
既存の描画処理にオフセットを加算して、縮尺の係数を掛けるだけで対応できた。
拡大して、特定の音域にスクロールして、マウスでハイライトすることで、採譜が行いやすくなった。
拡大縮小とスクロールを実装した。
これで、実装したい機能はほぼ実装できた。
Gumroadで配布しようと思っていたが、payoutの設定がうまくいかず別の手段を検討している。
その６で、低域をバンク分割することで、時間分解能を上げつつ、最小音階をC1まで下げることができた。
その際に、低域のバンクでは、オクターブあたりのビン数が高域の半分になるため、本来のCQTのビンの間隔の半分の間隔でスペクトルカーネルを作成することで、高域とビン間隔が同じになるように補間を行った。
これにより、
となり、オクターブ当たりの72グリッドに統一されることになる。
ついでに高域のバンクをさらに分割して、
のバンクを追加する。
3つのバンクで、高域のバンクほど周波数分解能が高くなる。
F2以下が低域バンクになっているため、ぼやけている。
これは、時間分解能とのトレードオフのため、補間を行っても本質的には改善しない。
時間分解能と周波数分解能のトレードオフをパラメータで設定できるようにすると良いかもしれない。
放物線補間を廃止して、スペクトルカーネルを本来のビン間隔より狭くすることで補間を行うようにした。
また、高域のバンクを追加して、高域の周波数分解能を上げた。
