前回は「実験材料の音源」である救急車のピーポー音をFFTしてピークとなる周波数の抽出を試みました。単純な最大値抽出だったので１周波数のみ検出。ピーポーで２波長、それにドップラー効果がかかっているので合計４波長が観察できるはずです。今回は「素朴なピーク検出関数」を自前で拵えて音源波形に適用してみたいと思います。
※「手習ひデジタル信号処理」投稿順 Indexはこちら
※Windows11上のScilab6.1.1およびScilab上のツールボックスを使用させていただいております。
前回作成のピーク検出関数のアップデート
前回作成のピーク検出関数はFFTかけて、その最大値のところのピーク値と周波数「だけ」を返す関数でした。今回は前回関数のFFT結果を踏まえた上でそのピーク値複数個を抽出したいので、前回作成関数の戻り値にFFT結果まで含めるようにし、その先に複数ピーク検出用の別関数を適用することにしました。
まずは前回のアップデート
// detect Peak FFT(2)
// Fs: sampling frequency [Hz]
// wav: time domain seq(nSample)
// opt: Window function 0:none, 1:hamming
function [pk, pfreq, xf, yDB] =detectPeak2FFT(Fs, wav, opt)
N = size(wav, '*');
regN = 2/N;
select opt
case 1 then
winf = window('hm',N);
regF = 2;
else
winf = ones(1,N);
regF = 1;
end
y = fft(wav .* winf);
xf = Fs*(0:(N/2))/N; //associated frequency vector
nf = size(xf, '*');
yDB = 20*log10(abs(y(1:nf)) * regN * regF);
pk = [];
pfreq = [];
[pk, wk] = max(yDB);
pfreq = xf(wk);
endfunction
前回結果のおさらい
前回のピーク１点のみ検出の結果のおさらいが以下に。上記のアップデートされたdetectPeak2FFTを使用。
Fs=22050;
y=ssAmbulance(Fs, Fs*6, 60, 0, 50);
[pk, pfreq] = detectPeak2FFT(Fs, y, 1)
clf;
plotFFT2(Fs, y, 1);
plot(pfreq, pk, 'x');
xtitle('peak detection(1)');
xgrid();
xlabel("freqency[Hz]");
ylabel("Mag[dB]");
結果は以下のように１点（×印のところ）のみ検出です。
もっとも素朴なピーク検出
もっとも素朴なピーク検出は微分して微分係数がプラスからマイナスに変じるところを捉えるってやつじゃないかと思います。学校で習うのは中学だった？あれ高校？忘却力の年寄には遥か太古の記憶です。元データの音源波形に対して、微分係数の正負変化点を×で示したものが以下に。
上記はデータ全体のうち着目部分のみなのですが、それにしてもピーク多過ぎ。これほど多いとどこを取り出したらよいのか分かりませんな。。。
素朴だけれども「も少し」絞り込んだピーク検出
上のグラフは局所的に計算しすぎて候補が多くなっているので「あるローカルな範囲内のチャンピオン」だけをピークとして選出するようにした関数が以下であります。NRに与える値により「ローカル」な範囲の広さを変更することが可能です。
// detect Peaks
// x: x vector
// y: y vector
// NR: range to search a local peak
function [pk, pfreq] =detectPeaks(x, y, NR)
N = size(y, '*');
if NR*2+1 > N then
error('ERROR: NR too large');
end
if NR < 2 then
error('ERROR: NR must be >2');
end
pk=[];
pfreq=[];
for idx = NR:N-NR+1
flag = %T;
for lidx = 1:NR-1
if y(idx) < y(idx-lidx) then
flag = %F;
break;
end
if y(idx) < y(idx+lidx) then
flag = %F;
break;
end
end
if flag then
pk($+1) = y(idx);
pfreq($+1)=x(idx);
end
end
endfunction
上記の関数を使って、再度ピーク検出を試みたものが以下に。
[pk3, pfreq3] = detectPeaks(xf, ydb, 11);
clf;
plotFFT2(Fs, y, 1);
plot(pfreq3, pk3, 'x');
xtitle('peak detection(3)');
xgrid();
xlabel("freqency[Hz]");
ylabel("Mag[dB]");
上記のグラフが以下に。
「裾」の方にはまだごちゃまんとピークが見えているのですが、このあたりの着目部分では「こんなもん」かも知れまへん。
なお、マジックナンバーである「１１」は、２から順番に適用していって「候補数」の減り方を調べてみた結果を見るという、アドホック？でヒューリスティック？な方法で決めました。でたとこ、やっつけ。
抽出結果はベクトルで返っているので、これをソートすれば上の方から結果を取り出せます。上位４つだけを取り出す操作はこんな感じ。
[pk3sort, sk] = gsort(pk3); pfreq3(sk(1:4))
上位４つの結果（周波数Hz）が以下に。
「ピー」音が約９１５Hzと約１００９Hz（ドップラー効果で２周波数に分かれて聞こえている）、「ポー」音が約８１０Hzと７３４Hzという結果が得られました。
まあ、ここでは「素朴な」アルゴリズムでなんとかなった感じ。でも、頭が平らなピークとかあるときには対応が必要だな。