前回「実験材料の音源」として救急車のピーポー音（サイレン）を作製。遠くから近づく様子が聞き取れて悦に入りました。さてその処理の第一歩として音源の周波数を求めたい、と。FFTかけてグラフ上で読み取れば一目瞭然なんでありますが、関数の処理結果で数値として取得したいです。とりあえずの関数作ってみたけれども機能不足デス。
※「手習ひデジタル信号処理」投稿順 Indexはこちら
※Windows11上のScilab6.1.1およびScilab上のツールボックスを使用させていただいております。
前回作成の「ピーポー音」関数のアップデート
早速のアップデートです。前回作成の版は音源としては使用できるのですが、結果を列ベクトルで返すようになってました。作成済の自前関数どもは行ベクトルに作用するようになってました。そのままだと処理の途中で行と列がマッチせず、エラーになる度に転置しないとならなくなりそうです。そこで結果は行ベクトルに統一することにいたしました。
// Sound Source: Ambulance
// Fs: sampling frequency [Hz]
// nSample: number of samples
// vs: sound source speed[km/h]
// v0: observer speed[km/h]
// dist: initial distance between source and observer [m]
// NOTE: Sound Level will be saturated under 1 [m]
function y=ssAmbulance(Fs, nSample, vs, v0, dist)
T06 = int(0.6 * Fs);
vsM = vs * 10^3 / 3600;
v0M = v0 * 10^3 / 3600;
AIR15=340.7;
dopplerAIR15=(AIR15-v0M)/(AIR15-vsM);
fqH = 960 * dopplerAIR15;
fqL = 770 * dopplerAIR15;
mag = 1.0;
cnt = 0;
flag = %T;
for idx = 1:nSample
timeNow = idx/Fs;
distNow = dist - (vsM-v0M)*timeNow;
if flag & (distNow < 0) then
dopplerAIR15=(AIR15+v0M)/(AIR15+vsM);
fqH = 960 * dopplerAIR15;
fqL = 770 * dopplerAIR15;
flag = %F;
end
if abs(distNow) < 1 then
distNow = 1;
end;
magNow = mag / abs(distNow);
cnt = cnt + 1;
if cnt > T06 then
cnt = 1 - T06;
end;
if cnt < 0 then
fq=fqH;
else
fq=fqL;
end;
yT(idx) = magNow * sin(2*%pi*fq*(idx/Fs));
end;
y = yT';
endfunction
上記関数により、以下条件で音データ y を作るシーケンスのが以下に（引数は前回とまったく同じ。）
-
- 前方50m にいる救急車がピーポーと鳴らしながら
- 時速60km/hの速度で
- 静止した観察者に向かって突進してきてその脇をすり抜ける
Fs=22050; y=ssAmbulance(Fs, Fs*6, 60, 0, 50);
参照のため、周波数777Hzの正弦波の音源波形 y777 を作るシーケンスが以下に。
[t777, y777]=genSin(Fs, Fs*6, 777, 0, 1);
上記で使っている自前関数 genSin() のソースも何度目かになりますが掲げておきます。
// Generate Sine wave
// Fs: sampling frequency [Hz]
// nSample: number of samples
// fq: frequency [Hz]
// phdeg: phase [deg]
// mag: manitude
function [t,y]=genSin(Fs, nSample, fq, phdeg, mag)
t = ((1:nSample) - 1) / Fs;
y = mag*sin(2*%pi*(fq*t+(phdeg/360)))
endfunction
ピーク周波数検出関数（第１版）
上記のようにして生成した「音」の周波数ピークを求めるための関数（第１版）が以下に。
// detect Peaks FFT
// Fs: sampling frequency [Hz]
// wav: time domain seq(nSample)
// opt: Window function 0:none, 1:hamming
function [pk, pfreq] =detectPeakFFT(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
上記の関数を使って、まずは 777Hzの音波形を処理してみたところが以下に。
[pk777, pfreq777] = detectPeakFFT(Fs, y777, 1)
処理結果が以下に。
777Hzが返ってきてますな。
つづいてピーポー音を食わせるシーケンスが以下に。
[pk, pfreq] = detectPeakFFT(Fs, y, 1)
その処理結果です。
「ピーポー音」のＦＦＴ結果のグラフが以下に。
「ピー」音と「ポー」音の２波長が、近づくときと遠ざかるときのドップラー効果で２波長に「分離」しているので、合計４波長あるわけね。グラフ上では上記のようにピークは一目瞭然だけれども、作成した関数では「ピー」音で近づいてくるときのピークのみ検出しているということみたい。
４種類全部を検出するためには「ピーク検出」に踏み込まないとダメなようだね。。。次回か。