R言語を用いたWavファイルの周波数解析

録音したWavファイルをR言語で周波数解析をおこないます。日本語がディレクトリにあるとエラー出ます。R言語をインストール後にRStudio Desktop（Free）をインストールします。ここでキャンパスで録音した10秒程度の虫の鳴き声（bugs.wav）を例にして解析を行います。
注意：日本語ディレクトリのエラーを避けるため2つのソフトをwindwosにインストールする際は、デフォルト指定の場所ではなくCドライブ直下などを選択ください。

使用 R version 4.3.1 (2023-07-6)
使用 RStudio 2023.06.0-421 (2023-07-6)

Sonic Visualiserによる音源の抽出
作業ディレクトリの変更
Wavファイルの読み込み
高速フーリエ変換（FFT）
スペクトログラム
バンドパスフィルタの設計と適用
バンドパスフィルタをかけた音源の出力

Sonic Visualiserによる音源の抽出

まずはSonic Visualiserをインストールしてから立ち上げてください。撮影した動画や音源をドラッグアンドドロップでファイルを読み込みます。次にSonic Visualiserのメニューから下記を実行してR言語で処理する音源を日本語を含まないディレクトリに保存してください。

使用 Sonic Visualiser (2023-07-6)
実行

File>Export Audio File...

作業ディレクトリの変更

以降はRStudioで操作します。まずは保存したWavファイルがあるディレクトリを作業ディレクトリに変更しましょう。RStudioのメニューから下記を実行して新しいプロジェクトフォルダを作成してください。

実行

Session＞Set Working Directory>Choose Directory...

Wavファイルの読み込み

インストールしたライブラリを使い読み込みます。今回はtuneRをインストールしてください。ライブラリは「Tools＞Install Packages...」で、PackagesでtuneRを選んでインストールできます。インストール後、下記のコマンドでtuneRを登録して、audioに保存したwavファイル(bugs.wav)を入れます。それぞれ文末で実行するためにエンターキーを押してください。

コンパイル実行例

library(tuneR)
audio=readWave("bugs.wav")

Wavファイルの内容確認は以下の通りstr(audio)と入力して実行します。

プログラム実行例

str(audio)

Formal class 'Wave' [package "tuneR"] with 6 slots
..@ left : num [1:481280] 0 0 0 0 0 0 0 0 0 0 ...
..@ right : num(0)
..@ stereo : logi FALSE
..@ samp.rate: int 44100
..@ bit : int 32
..@ pcm : logi FALSE

Wavファイルの情報が続けて表示されます。ここではleftに音源が入っており、rightに入っていないのでモノラル音源だと分かります。また、samp.rateからサンプリング周波数が44100 Hzです。次に以下のコマンドを実行して10秒の音源から4秒から6秒の波形をグラフで表示します。

4秒から6秒の波形をグラフで表示

f_s <-audio@samp.rate
audioL=audio@left
s=f_s*4
t=f_s*6
plot((s:t)/f_s,audioL[s:t],type="l",xlab="time [s]",ylab="Signal [V]")

図1にグラフが表示されます。背景音もありはっきりした特徴はわかりにくい波形です。

図1：4秒から6秒までの波形

上記のグラフ出力はplotのコマンド前後に下記のようにコマンドを追加することで可能です。

プログラム実行例

png("plot.png",width=500,height=500,pointsize =18)
plot((s:t)/f_s,audioL[s:t],type="l",xlab="time [s]",ylab="Signal [V]")
dev.off()

ページの先頭へ

高速フーリエ変換（FFT）

図1でみた時系列データを高速フーリエ変換により特徴的な周波数を調べます。下記のコマンドを実行していきます。

時系列データの高速フーリエ変換（FFT）

amp.AudioL <- fft(audioL)
amp.AudioL.cmp <-Mod(amp.AudioL)
plot(amp.AudioL.cmp,type="l")

上からコマンドを見ていくと、audioLをFFT関数で結果をamp.AudioLに入力します。Mod関数で複素振幅から振幅の大きさを計算します。まずはplotで図を表示します。縦横の軸は特に指定していませんが、横軸の真ん中で左右対称になっていることがわかります。

図2：bugs.wavのFFT結果

ナイキスト周波数で折り返した左半分のグラフを画像でフォルダへ出力します。

ナイキスト周波数で折り返した左半分のグラフ

f <- seq(from=0, to=f_s/2, length=length(audioL)/2+1)
AudioL.ampspec <- amp.AudioL.cmp[1:(length(audioL)/2+1)]
png("plot.png",width=500,height=500,pointsize =18)
plot(f, AudioL.ampspec, type="l", xlab="Frequency (Hz)", ylab="Amplitude Spectrum")
dev.off()

図3：bugs.wavのスペクトル分布

今回は横軸に周波数を表示しています。スペクトルをみると、0 Hzの低周波領域に強い背景音の成分と3000 Hzから6000 Hzの間に特徴的な2つのピークがあります。

ページの先頭へ

スペクトログラム

フーリエ変換をすると残念ながら時間の情報がなくなります。そこで短い時間領域をフーリエ変換して、大まかな時間毎に周波数変化をみていきます。この操作は短時間フーリエ変換と言い、得られるグラフはスペクトログラムと呼ばれます。ドラマや映画の科学捜査の一場面で写る声紋分析などにつかわれます。時間領域を窓サイズと呼び、またオーバーラップさせながら描画します。サイズやオーバーラップは波形に合わせて調整します。ここから、スペクトログラムを使うのでsignalとカラー画像で出力するためにoceの2つのライブラリを使うのでtuneRと同様に下記コマンドを実行する前にインストールしてください。signalはRStudioの再起動が必要になります。

スペクトログラムの実行

library(signal, warn.conflicts = F, quietly = T)
spec=specgram(x=audioL,n=1024,Fs=f_s, window=256, overlap=128)
P = abs(spec$S)
P = P/max(P)
P = 10*log10(P)
t = spec$t

specgramの中で窓サイズとオーバーラップを指定しています。その他の詳細についてはここを一読ください。以下は画像でスペクトログラムをフォルダへ画像出力します。使われるPCの性能に依りますが少し時間が必要です。

スペクトログラムの出力

library(oce)
png("plot_specgram.png",width=500,height=500,pointsize =18)
imagep(x = t,y = spec$f,z = t(P),col = oce.colorsViridis,ylab = 'Frequency [Hz]',xlab = 'Time [s]',drawPalette = T,decimate = F)
dev.off()

図4：bugs.wavのスペクトログラム

図3と同様に4000 Hzと5000 Hz付近で特徴的なピークがあることがわかります。しかし得られたスペクトログラムでみると、5000 Hz付近で間欠的に4000 Hz付近で連続的に2匹の虫が鳴いているようです。もう一度音源を聴いてみてもやはり背景音が大きいためよく分かりません。

ページの先頭へ

バンドパスフィルタの設計と適用

そこで、3000 Hzから6000 Hzの間の成分だけを取り出すバンドパスフィルタの設計と適用をおこないます。フィルタはfir1関数を用いていて、詳細はここを参考にしてください。

バンドパスフィルタの設計

fil_N <- 512
fn <- f_s/2
fc <- c(3000,6000)
fc_norm <- fc/fn
fir_filter <- fir1(fil_N, fc_norm, type="pass")
freqz(fir_filter, Fs=f_s)

最後のコマンドfreqzで設計したフィルタを表示しています。fil_Nとfcをいじるとフィルタの特性が変わります。以下にフィルタを音源に適用してフォルダへ画像を出力します。

バンドパスフィルタの適用

audio_bpfil <- filtfilt(fir_filter, audioL)
png("plot_b_signal.png",width=500,height=500,pointsize =18)
plot((s:t)/f_s,audio_bpfil[s:t],type="l",xlab="time [s]",ylab="Signal [V]")
dev.off()

図5：3000 Hzから6000 Hzの間の成分だけを取り出したbugs.wavの信号

図1とくらべると埋もれていた波形がみえてきました。スペクトログラムでも下記のコマンドを実行することでバンドパスフィルタの適用を確認できます。

バンドパスフィルタの適用

spec=specgram(x=audio_bpfil,n=1024,Fs=f_s, window=256, overlap=128)
P = abs(spec$S)
P = P/max(P)
P = 10*log10(P)
t = spec$t
#output graph in png file
png("plot_b_specgram.png",width=500,height=500,pointsize =18)
imagep(x = t,y = spec$f,z = t(P),col = oce.colorsViridis,ylab = 'Frequency [Hz]',xlab = 'Time [s]',drawPalette = T,decimate = F)
dev.off()