MATLAB mobileを用いて電車の速度を測る

はじめに

　年１回しか投稿しないという、このQiita Advent Calender向けのMATLAB記事、今年も出番が回ってまいりました。
　昨年までは、MATLAB-Excel連携のハウツーを書いてきたのですが、どうもAIに押されて芳しくない状況になってきていますので、読み物的な内容を目指そうと思います。

１．よくあるアプローチ

　MATLAB mobileで速度を求めるアプローチとして、誰もが考えるのはGPSデータ（位置データ）を取得する手法かと思います。下図にもありますように何も考えないで（失礼）、速度を見ることができます。

２．今回のアプローチ

　かつて、鉄道の線路は短いレールを継ぎ合わせて形成されていました。一般的には25m長のものが多いようです。それで、乗車しているとガタンゴトンという音がして、大まかな速度は把握できたものです。例えば時速80kmで走行している場合、「ガタンゴトン」の周期は、およそ1.1秒となります。
しかし、今はロングレールという接合されたレールが使われていることが多く、ガタンゴトンという音は聞けなくなっています。
ですが、諸般の事情により、継ぎ目が残してある部分があります。こういうところを通過するときには、ガタンゴトンという音がします。
ただし、（１）線路脇で観測した場合と、（２）電車内で観測した場合では聞こえ方が異なります。
（１）線路脇で観測した場合
　　この場合は、全ての車両の車輪が、１カ所の継ぎ目において音を発生するので、「ガタン、ガタンゴトン、（略）、ゴトン」というように聞こえます。
（注）車両が標準的な長さ（20m）の場合、時速80kmで通り過ぎたとすると、「ガタンゴトン」の周期は0.9秒となります。つまり電車内で観測した場合に較べて周期は短くなります。
（２）電車内で観測した場合
　　この場合には、自分が乗っている車両の車輪が継ぎ目を通過するときの音のみとなります。
　「ガタン」もしくは「ゴトンガタン」の単発１回のみ。
ここで、今回は、この「ガタン」１回のみから速度を求めてみようというものです。
電車の車体には前後に１つずつの台車がついており、台車には前後に２軸の車輪がついています。１つの台車の２つの車輪の距離（軸距離、クルマで言うところのホイールベース）は、首都圏の一般的な車両においては2100mmとなっています（京王線はちょっと違うらしい）。
この2100mmの軸距離で「ガタン」という音がするわけですが、「ガ」と「タン」の間の時間は、時速80kmだと、94.5msecとなり、人の感覚だと短すぎてピンと来ない領域になります。
そこで、今回、これをMATLAB Mobileで求めて、時速に換算してみようということです。

３．とりあえずの音源解析

　線路に継ぎ目があるところに行って、録音してきました。ここは鉄橋があるので継ぎ目があります（線路の伸び縮みの関係？　ちょっとそこはわからない。ちなみにトンネル内も継ぎ目がありますね）。
録ってきた音をMATLABに読み込ませて、波形を表示します。

[A,fs]=audioread("Sound01.mp4");
figure(1);
t=(1:length(f))/44100;
plot(t,f(:,1));
xlabel("time[sec]");
ylabel("Amplitude"); 

　「ガタン」「ゴトン」の部分は、わりとわかりやすそうです。車両数も６両あることがわかります。振幅がだんだん小さくなっているのは、自動ゲイン制御でしょうか。

４．波形の抽出

　さて、この波形から「継ぎ目音」の部分を抽出してみます。
１両目の先頭部分の音（上の図の①の左部分）を拡大した図を上段に、そのスペクトログラムを中段に示します。
低い周波数においては、「継ぎ目音」と関係なく、音がしていることがわかります。そこで、ハイパスフィルタを用いて低い音をカットしたのが下段の波形となります。カットオフ周波数は4kHzとしました。

次に、この波形の包絡線を計算します。ここではsignal proccesing toolboxのenvelope関数を使ってみました。
下の図の波形において、上段がハイパスフィルタをかけない場合、下段がハイパスフィルタをかけた場合です。ハイパスフィルタをかけたほうが「継ぎ目音」の抽出がしやすくなっていることがわかります。

「継ぎ目音」のタイミングを抽出します。
ここでは、signal processing toolboxのfindpeaks関数を使っています。findpeaks関数は、パラメータを指定しないと、あらゆるピーク部分を拾ってしまうので、狙った部分が取れないかと思います。
今回はピークのプロミネンスを指定して、該当箇所を絞っています。

このようにして得た、継ぎ目音のタイミングから速度を求めます。
下の図のように、全車両分の波形データにハイパスフィルタをかけ、包絡線を計算し、継ぎ目音部分の推定を行ないます。
継ぎ目音のタイミングの差分を取り、そのうちの一番時間が短いものを探します。ここは、台車の２つの車輪間ということになりますので（となりの台車の車輪までの距離のほうが長いため）。

そして、台車の２つの車輪間（2.1m)と音の時間差を元に、時速を求めグラフに表示します。
ここでは、おおよそ時速40kmという結果が算出できています。

５．MATLAB mobileを用いて音を取りこみ測定する

MATLAB mobileを用いてスマホのマイクで音を取得するサンプルについては、MathWorksのここのページが参考になります。

モバイル デバイスからのオーディオ データの取得と解析

上記のページを参考にスクリプトを書いてみました。

% モバイルデバイスの作成
mobileDevObject = mobiledev;
% マイクロホンの許可
mobileDevObject.MicrophoneEnabled = 1;

% バックグランド音のレベルを把握する
mobileDevObject.logging = 1;
disp('Sampling background noise...')
pause(2)

% データ取りこみ停止
mobileDevObject.logging = 0;

% モバイルからオーディオデータを吸い上げる
audioData = readAudio(mobileDevObject);

% フィルタリング
fc=4000;
fs = mobileDevObject.Microphone.SampleRate; %サンプリング周波数
[b1,a1]=butter(5,fc/(fs/2),'high'); %ハイパスフィルタ
A2=filtfilt(b1,a1,audioData(:,1));

% 包絡線を計算
fl=500;
[YA,~]=envelope(A2,fl,'rms');

disp('Maximum sound of the background noise: ')
threshold = max(YA, [], "all")

% 継ぎ目音の検出
disp('Rail Sound detecting for a few seconds.')
% データ取りこみの開始
mobileDevObject.logging = 1;
tic
disp('"Recording MATLAB Mobile Audio"')
startTime = 0;
totalAudio = [];

% realaudioを200msecずつ録り、閾値の1.5倍のデータがある場合は
% startTimeをゼロでなくする
% 取りこみ時間は最大で30秒

while toc < 30 && startTime == 0
    pause(.2); % 200msec
    audioData = readAudio(mobileDevObject);
    A2=filtfilt(b1,a1,audioData(:,1));
    [YA,~]=envelope(A2,fl,'rms');
    % 閾値を越えたデータが来たとき
    if max(YA) > threshold * 1.5
        threshold2 = max(YA);
        startTime = toc;
        totalAudio = audioData;
    end
end

% レール音が終了するまで録りつづける
if startTime ~= 0
    numberOfIntervalsStopped = 0;
    while numberOfIntervalsStopped < 2 && toc < 10
        pause(.2)
        audioData = readAudio(mobileDevObject);
        A2=filtfilt(b1,a1,audioData(:,1));
        [YA,~]=envelope(A2,fl,'rms');
        if max(YA) < threshold * 1.5
            numberOfIntervalsStopped = numberOfIntervalsStopped + 1;
        else
            numberOfIntervalsStopped = 0;
        end
        totalAudio = vertcat(totalAudio,audioData);
    end
end
% 取りこみ終了
mobileDevObject.logging = 0;

endTime = toc;

leftAudio = totalAudio(:,1);
n = numel(leftAudio);
if n == 0 % 音が撮れなかった場合
    disp(' ')
    disp('No audio data recorded. Try to run the script again.')
    clear mobileDevObject
    return
end

%%%% 音が入ってきた場合の速度解析ルーチン
disp('Maximum sound level[dB]: ')
disp(20*log10(abs(threshold2/threshold)));

figure(1);
plot(leftAudio)
title('Sound wave');
timeElapsed = endTime - startTime
ticks = 0:floor(timeElapsed);
sampleTicks = ticks * n/timeElapsed;
xticks(sampleTicks)
xticklabels(ticks)
xlabel('Time(s)')
ylabel('Amplitude')

%%%%%%%%%%%
fc=4000;
[b1,a1]=butter(5,fc/(fs/2),'high');
A2=filtfilt(b1,a1,leftAudio);
leftlen=length(leftAudio);
tt=(1:leftlen)/fs;
A3=A2;
figure(3);
subplot(3,1,1);
plot(tt,A3);

fl=500;
[YA,~]=envelope(A3,fl,'rms');
subplot(3,1,2);
plot(tt,YA);

[YA3,~]=envelope(A3,fl,'rms');
subplot(3,1,3);
plot(tt,YA3);
MPP=0.003;
findpeaks(YA3,tt,'MinPeakProminence',MPP);%,'Annotate','extents');
[pk,lk]=findpeaks(YA3,tt,'MinPeakProminence',MPP,'Annotate','extents');

%%% 速度計算
Velocity=2.1*3600/1000/(lk(2)-lk(1))
fc=4000;
[b1,a1]=butter(5,fc/(fs/2),'high');
A2=filtfilt(b1,a1,A(:,1));
A3=A2;
figure(4)
subplot(2,1,1);
plot(t,A3);
xlabel("time[sec]");
ylabel("Amplitude");

[YA3,~]=envelope(A3,fl,'rms');
subplot(2,1,2);
plot(t,YA3);
% plot(YA3);
findpeaks(YA3,t,'MinPeakProminence',MPP);%,'Annotate','extents');
[pk,lk]=findpeaks(YA3,t,'MinPeakProminence',MPP,'Annotate','extents');
DF=diff(lk);
DFlen=length(DF);
Odd=rem(DFlen,2);
if Odd
    DP=1:2:DFlen;
else
    DP=2:2:DFlen;
end
DFtime=DF(DP);
Velocity=round(2.1*3600/1000./DFtime);
text(lk(DP),DP*0,num2str(Velocity'),'FontWeight','Bold','color','red');
xlabel("time[sec]");
ylabel("Amplitude");

clear mobileDevObject

6．電車の中で測ってみる

やってみました。
電車の中でおしゃべりしている人は少なく、それはいいのですが、アナウンスが流れている時間が相当長いです。最近は英語アナウンスもありますし。
なので、「うまくすれば速度が測れる」くらいの感触です。
うまくいった例を下記に示します。

時速83kmです。

現場からは以上です。