手習ひデジタル信号処理(131) Scilab、素朴なピーク検出、その2

Joseph Halfmoon

前回、自前の「素朴なピーク検出関数」を作製、FFT結果に適用し所望のドップラー効果のピーク周波数位置を求められました。しかし、ありがちな「ピークが台地(メサ)状だったらどうよ」問題に対する挙動を確かめてません。お惚け老人は自分で作っておきながら既にロジックは忘却力の彼方。実際に確かめてみるしかないっと。どうなのよ。

※「手習ひデジタル信号処理」投稿順 Indexはこちら

※Windows11上のScilab6.1.1およびScilab上のツールボックスを使用させていただいております。

「ピークが台地(メサ)状」

大陸では、巨大なmasaが聳え立っていたりしますが、日本じゃ荒船山くらい?MESAと言えば頭が平らな山というか、台地といった地形です。ここで「素朴なピーク検出アルゴリズムの天敵」として登場するのが下のような波形ね。mesa

頭が平らっす。素朴なピーク検出アルゴリズムは高校(中学校?)でならう「微分係数が+から‐への変化点をピークとする」を基本にしているので、単純に適用するとこの手の「メサ状地形」を見逃すことになります。台地とは言えピークはピークであるハズなのだけれど。

まずはテスト用波形の生成関数作成

まずは上記のようなメサ状のテスト波形を作らないとなりません。元波形は「いつもの」正弦波生成関数で生成するものとします。

Fs=22050;
nSample=22050;
freq=777;
[t0, y0]=genSin(Fs, 1000, freq, 0, 1);

上記で生成した正弦波 y0 の頭のところをつぶして台地状にする自前関数を作製してみました。こんな感じ。

// upper limit
// x: input vector
// lev: limit, input value ceilled at this level
function xlim = upperLimit(x, lim)
    xlim = x;
    N = size(x, 2);
    limV = linspace(lim, lim, N);
    xlim(x > limV) = lim
endfunction

limitを超えるポイントをサチらせるところは for を使わずにScilab風に書けて良かったですが、入出力とも 1 x N サイズの「行ベクトル」形式なのはいかがなものか。最初に作った genSin関数が「行ベクトル」形式なのでそれに合わせているのだけれども、Scilab的(数学的)には列ベクトルの方が良かった気がします。後の祭りだけれども。

上記の関数を使って、正弦波の頭を0.8で平らかにするのが以下に。

ylim = upperLimit(y0, 0.8);
前回作成した「素朴なピーク検出アルゴリズム」を適用してみた

上で作った波形に前回のピーク検出アルゴリズムを提供してみるシーケンスが以下に。後ろの方は結果のプロット用ね。11というのは前回FFTの時に使ったマジックナンバーっす。ぶっちゃけ処理対象によって「出し入れする」部分。

[pk0, pt0] = detectPeaks(t0, y0, 11);
[pklim, plt0] = detectPeaks(t0, ylim, 11);
clf;
plot(t0', [y0', ylim']);
plot(pt0, pk0, 'x');
plot(plt0, pklim, 'x');

処理結果のプロットが以下に。青がオリジナルの正弦波、緑が頭がメサな波形、×印が検出されたピークです。testWaveAndPeaks

上記じゃよくわからないので頭の部分の拡大が以下に。detectedPeaks

×印を見てみると、危惧されたメサ状の「フラットトップ」な部分もピークとして認識してはいるみたいです。しかし、時間的に先頭の赤丸つけたところはピークとして認識してません。これは例の「マジックナンバー」、ピークを探すローカル範囲の大きさのせいです。系列の冒頭および末尾のマジックナンバー数より少ないサンプリングポイントを無視してしまうためです。これは「マジックナンバー」使っている限りいたしかたないのう。。。

これからするとメサ状のピークは検出できる感じ。けれども「踊り場、肩」みたいな一端フラットになってまた上昇したりするポイントは、「マジックナンバー」次第で、純粋なピークではないのだけれども検出してしまうっと。まあ、分かった上で使う分には使えるようだが、忘却力だからな、忘れそうだな。。。

手習ひデジタル信号処理(130) Scilab、素朴なピーク検出、その1 へ戻る

手習ひデジタル信号処理(132) プロは1行でピーク検出。Scilabでなぞってみた へ進む