Excelで簡易的にWAVファイルから周波数分析をしてみよう
こんにちは。すごい改善のAitoです。
今日は少し変わったネタで書いていきます。それも「ExcelでWAVファイルから周波数分析を行う」です！一体何を言っているんだ？という方も、これやったことあるという方も様々な反応が出ると思います。理系で学んでいると専攻によってはどこかで音声処理を学ぶ機会があります。ここでは最も身近な分析ソフトであるExcelを使って解析をしてみましょう。ロジックの勉強にもなって一石二鳥です。少し重ための内容になりますが、是非お付き合いください！
では始めていきましょう！
まずは、事前知識から。
さて、音とは何でしょうか。音は空気の振動で、疎密波と呼ばれる波です。例えばこんな感じのものです。
実際に私たちが聞く音は、連続した波形になります。いわゆるアナログデータというやつですね。一方、コンピュータで扱う値はごく短い時間で区切って音の強弱を数値で表現した離散値を取ります。こちらはいわゆるデジタルデータと呼ばれるものですね。上記の図はデジタルデータをAudacityというソフトで可視化したものになります。
では、この音を分析していくにはどうしましょう？巷ではよく「データ分析！」と叫ばれることが増えましたが、ではどうやって分析するの？となると途端に「わからない！」となる方も多いのではないでしょうか。一度やり方を見てしまえば、自分でもできるようになると思いますので、今回はこんな風にやるんだという一例として見てみてください。話がそれました、分析手法についてですが、今回は「(離散)フーリエ変換」を用います。フーリエ変換とは簡単に言うとある音(の波形)の中にどのような周波数が含まれているかを調べる手法です。入力信号(いわゆるこちらを表示したものが波形になります)の種類で分析手法が変化します。アナログで周期的な信号の場合はフーリエ級数展開を用い、アナログで周期的な信号の場合はフーリエ変換を用い、デジタル信号の場合は離散フーリエ変換を用います。今回はデジタル値を扱いますので、離散フーリエ変換ですね。
離散フーリエ変換について少し説明しておきます。波を表す関数は
と表すことができます。関数は三角関数の和の形で表すことができると言えます(オイラーの公式より、三角関数に変換可能)。係数は各周波数における周波数成分の大きさを表しているので、2番目の式に関してkの関数F(k)とみなすと、
となります。デジタル値は離散値をとるので、積分を総和に置き換えてあげると、
と表せます。これで各に対応する周波数成分の割合を特定することができます。
今回プログラムで実装するにあたって、この離散フーリエ変換を限定的な場面で高速に実行できる高速フーリエ変換を用いますが、原理は同じだと考えてもらって大丈夫です。この実装はこちらの記事を参考にして使わせていただきます。
【FFT】
http://blog.livedoor.jp/sce_info3-craft/archives/7774799.html
次はコンピュータにおいて、音の取り扱いについて解説します。マイクに入力された音はごく短い時間で音量が正規化されて振幅が保存されます。このデータを保存しておくのがWAVファイルになります。WAVファイルは音の波形が何も加工されずに直接書かれている形のファイルです。最近はあまり見なくなってしまった気がしないでもないですが、一番音質が良い形式の一つとなります。(いわゆる無劣化保存ができます)
Excelでこのファイルを読み込むにはVBAを使わないといけません。こちらも先人がありがたいことに読み込みのための記事を書いてくださっていたので、少し手を加えて使用します。
【WAVファイル読み込み】
http://otktake.blogspot.com/2015/01/excelwave.html
さて、基礎的なところは解説したので、次は実際に読み込んで解析するところに入っていきましょう。
まずはExcelを起動し、Alt + F11を押してVBEを起動します。起動したら上の挿入メニューから標準モジュールとクラスモジュールを追加し、次の図のように名前を変更します。
ここまで出来たら次はWAVファイルの読み込みから作っていきます。今回解析する曲は「カエルの歌」にでもしておきます(笑)
ここからダウンロードしてください。この音源を使用します。
先ほどの記事を参考に、書きます。VarModuleに次のコードを書きます。
'共通変数
'**************************************************************WAVEファイル***************************************************************************************
Public Type RIFF_CHUNK
ChunkID As String * 4 '"RIFF"(0x52494646)で固定
ChunkSize As Long 'チャンクサイズ。ファイル全体のサイズからファイル全体サイズからRIFFとWAVEのバイト数(8Byte)を引いた数。この情報をもとにWavファイルのファイルサイズを算出できる。
Format As String * 4 '"WAVE"(0x57415645)で固定
End Type
Public Type fmt_chunk
SubchunkID As String * 4 '"fmt"(0x666D7420)で固定
SubchunkSize As Long 'fmtチャンクのバイト数。リニアPCMならば16(0x10000000)。その他は16+拡張パラメータ
AudioFormat As Integer '音声フォーマット。非圧縮のリニアPCMフォーマットは1(0x0100)。
NumChannels As Integer 'チャンネル数。モノラルは1、ステレオは2
SamplingRate As Long 'サンプリング周波数(Hz)。44.1kHzの場合なら(0x44AC0000)
ByteRate As Long '一秒あたりのバイト数の平均。サンプリング周波数*ブロックサイズで求める。44.1kHz、16bit、モノラルならば44100x2x1=88,200(0x88580100)
BlockAlign As Integer 'ブロックサイズ。チャンネル数 * 1サンプルあたりのビット数 / 8で求める。モノラル16bitなら1*16bit = 16bit = 2byte(0x0200)、ステレオ16bitなら4(0x0400)
BitsPerSample As Integer 'ビット/サンプル。1サンプルに必要なビット数。8ビットの場合は8(0x0800)、16ビットの場合は16(0x1000)など。
End Type
Public Type DATA_CHUNK
SubchunkID As String * 4 '"data"(0x64617461)で固定
SubchunkSize As Long '波形データのバイト数。
End Type
'******************************************************************************************************************************************************************
そして、WAVクラスモジュールに次のコードを書き込みます。
'**************************************************再定義**************************************************************************
Private riffChunk As RIFF_CHUNK
Private fmtChunk As fmt_chunk
Private dataChunk As DATA_CHUNK
Private data() As Integer
'**********************************************************************************************************************************
'シートに音階と周波数の対応表がある
Private Pitch_Hz(87) As Double
Private Pitch_Code(87) As String
Const PI = 3.14159265358979
Private Sampling As Long
Private Peek_Hz_Before As Long
'WAVファイルを読み込む(44.1kHz/16bit 44.1kHz/8bit 両方に対応)
Public Function Read_WaveFile(ByVal WaveFilePath As String) As Integer()
Open WaveFilePath For Binary As #1 'WAVファイルをバイナリ形式で開く
Get #1, , riffChunk 'RIFF識別子の読み込み
Get #1, , fmtChunk 'fmt識別子の読み込み
Get #1, , dataChunk '波形データに関する情報の読み込み
ReDim data(dataChunk.SubchunkSize / 2) As Integer 'dataの配列数をファイル形式に合わせて再定義
Get #1, , data '波形データの読み込み
Close #1
Sampling = fmtChunk.SamplingRate
Read_WaveFile = data '波形データを返す
End Function
これで読み込みができるようになりました。試しに読み込んで表示するプログラムを書いてみましょうか。
MainModuleに次のように書きます。
Sub ReadTest()
On Error Resume Next
Dim CWAV As New WAV 'WAVクラスモジュールを使う
Const FilePath As String = "D:\Song.wav" 'ファイルパス
Dim data() As Integer: data = CWAV.Read_WaveFile(FilePath) 'データの読み込み
Dim i As Long '表示用
For i = 0 To 4410 '0.1秒分表示(サンプリング周波数が44100Hzなので)
Cells(i, 1).Value = data(i) 'とりあえずA列に書き出す
Next
End Sub
すると、A列に振幅が読み込まれます。読み込まれた範囲をグラフ化してみましょう。
ちゃんと波形が読み込まれてますね！確認が取れましたので次は実際にフーリエ変換して周波数を分析してみましょう。
'ピークの周波数を返す
Public Function Get_Peek_Hz(ByRef FFTData() As Double, ByVal N As Long) As Double
'ピークピッキングによってピークを求める
Dim peaks_val() As Double
Dim peaks_index() As Long
Dim max_index As Long
Dim Counter As Long: Counter = 0
Dim i As Long: For i = 2 To UBound(FFTData)
If (FFTData(i - 1) - FFTData(i - 2)) >= 0 And (FFTData(i) - FFTData(i - 1)) < 0 Then
ReDim Preserve peaks_val(Counter) As Double
ReDim Preserve peaks_index(Counter) As Long
peaks_val(Counter) = FFTData(i - 1)
peaks_index(Counter) = (i - 1)
Counter = Counter + 1
End If
Next i
'最大値を求める
max_index = Get_MaxIndex(peaks_val)
max_index = peaks_index(max_index)
Dim ret As Long
If Peek_Hz_Before <> 0 Then
ret = (((Sampling / N) * (max_index)) + Peek_Hz_Before) / 2
Else
ret = (Sampling / N) * (max_index)
End If
Peek_Hz_Before = ret
Get_Peek_Hz = ret
End Function
'配列の最大値を含んでいる要素数を求める
Public Function Get_MaxIndex(ByRef TargetArray() As Double) As Long
Dim ret As Long: ret = -1
Dim tmp As Double: tmp = -1
Dim i As Long: For i = 0 To UBound(TargetArray)
If TargetArray(i) > tmp Then
tmp = TargetArray(i)
ret = i
End If
Next i
Get_MaxIndex = ret
End Function
'FFTをかけて音声の基本周波数を取得する
Public Function Get_Pitch_Hz(ByVal WaveFileName As String) As Long
Dim Num As Long: Num = 4096 'FFTのサイズ
Dim SoundData() As Integer: SoundData = Read_WaveFile(ThisWorkbook.Path + "\tmp\" + WaveFileName + ".wav")
Dim InputData() As Double: ReDim InputData(Num) As Double
Dim N, k As Integer
Dim x_real() As Double, x_imag() As Double
ReDim x_real(Num) As Double
ReDim x_imag(Num) As Double
'波形
For N = 0 To Num - 1
InputData(N) = CDbl(SoundData(N))
Next N
Call FFT(InputData, x_real, x_imag, Num)
'出力
Dim Out() As Double: ReDim Out((Num / 2) - 1) As Double
For k = 0 To (Num / 2) - 1
Out(k) = Sqr((x_real(k)) ^ 2 + (x_imag(k)) ^ 2)
Next k
Get_Pitch_Hz = Get_Peek_Hz(Out, Num)
End Function
Public Function FFT(x() As Double, ByRef y_re() As Double, ByRef y_im() As Double, N As Long) As Boolean
'引数
'x 入力データ
'y_re 戻り値実部
'y_im 戻り値虚部
'N FFTサイズ
'戻り値
'TRUE FFT成功
'FALSE データ数が2のべき乗でないので失敗
'変数
Dim PI As Double 'π
Dim StageCount As Integer 'ステージ数
Dim BlockCount As Long 'ブロック数
Dim NodeCount As Long 'ノード数
Dim stage As Integer 'ステージ番号
Dim block As Long 'ブロック番号
Dim node As Long 'ノード番号
Dim n1 As Long '計算用
Dim n2 As Long '計算用
Dim i As Long '計算用
Dim r As Double '計算用
Dim Index() As Long 'インデックス並べ替え用
Dim re1 As Double '計算用
Dim im1 As Double '計算用
Dim re2 As Double '計算用
Dim im2 As Double '計算用
'FFTサイズ確認
If CLng(N And (N - 1)) <> 0 Then
FFT = False
Exit Function '2のべき乗でなければ終わり
End If
'計算準備
PI = Math.Atn(1) * 4 'π
StageCount = CInt(Round(Log(N) / Log(2), 0)) 'ステージ数
y_re = x '実部
ReDim y_im(N) '虚部
'バタフライ演算の繰り返しループ
For stage = 0 To StageCount - 1
BlockCount = 2 ^ stage 'ブロック数
NodeCount = N / BlockCount 'ノード数
r = 2 * PI / N * BlockCount '定数
For block = 0 To BlockCount - 1
For node = 0 To NodeCount / 2 - 1
'バタフライ演算
n1 = node + NodeCount * block
n2 = n1 + NodeCount / 2
re1 = y_re(n1)
im1 = y_im(n1)
re2 = y_re(n2)
im2 = y_im(n2)
y_re(n1) = re1 + re2
y_im(n1) = im1 + im2
y_re(n2) = (re1 - re2) * Cos(r * node) + (im1 - im2) * Sin(r * node)
y_im(n2) = -(re1 - re2) * Sin(r * node) + (im1 - im2) * Cos(r * node)
Next node
Next block
Next stage
'並び替え処理
ReDim Index(N)
For stage = 0 To StageCount - 1
For i = 0 To 2 ^ stage - 1
Index(2 ^ stage + i) = Index(i) + 2 ^ (StageCount - stage - 1)
Next i
Next stage
For i = 0 To N - 1
If Index(i) > i Then
re1 = y_re(Index(i))
im1 = y_im(Index(i))
y_re(Index(i)) = y_re(i)
y_im(Index(i)) = y_im(i)
y_re(i) = re1
y_im(i) = im1
End If
Next i
FFT = True
End Function
Public Sub Class_Initialize()
'シートから読み込む
Dim i As Long: For i = 1 To Worksheets(1).Cells(Rows.Count, 1).End(xlUp).Row
Pitch_Hz(i - 1) = Worksheets("Sheet1").Cells(i, 1).Value
Pitch_Code(i - 1) = Worksheets("Sheet1").Cells(i, 2).Value
Next i
End Sub
'ピッチコードの取得
Public Function Get_PitchCode(ByVal Hz As Long) As String
Dim ret As String
Dim Difference(87) As Double
Dim i As Long: For i = 0 To UBound(Pitch_Hz)
Difference(i) = Abs(Hz - Pitch_Hz(i))
Next i
Dim MinimumIndex As Integer: MinimumIndex = Min_Index(Difference)
'いじった
' If Abs(Pitch_Hz(MinimumIndex) - Hz) < 3 Then
' ret = Pitch_Code(MinimumIndex)
' Else
' ret = "None"
' End If
If Abs(Pitch_Hz(MinimumIndex) - Hz) < 10 Then
ret = Pitch_Code(MinimumIndex)
Else
ret = "None"
End If
Get_PitchCode = ret
End Function
'配列の最小値をとるインデックスを返す。
Private Function Min_Index(ByRef TargetArray() As Double) As Long
Dim Index As Long: Index = 0
Dim min As Double: min = TargetArray(0)
Dim i As Long: For i = 1 To UBound(TargetArray)
If TargetArray(i) < min Then
min = TargetArray(i)
Index = i
End If
Next i
Min_Index = Index
End Function
このようにFFTとそれを扱う周辺のプログラムを書きます。これだけだと何をやっているかさっぱりだと思うので、内容の解説をしますね。フーリエ変換を行うと「Get_Pitch_Hz」プロシージャにFFTのサイズが4096と指定してあります。このサイズが重要でして、周波数分解能(どれくらい細かく周波数を調べられるか)が変化します。式は次のように与えられます。
つまりFFTのサイズが置きくなれば分解能は上がりますが、計算に時間がかかるようになります。今回はサンプリング周波数が44100なのでだいたい10.8Hzくらいの分解能です。ちょっと物足りない気がしますが、今回はゆっくりした曲(カエルの歌(笑))なので、まあ妥協してこれくらいでいいでしょう。さて、分析すると10.8Hz,21.6Hz,…のようにそれぞれの周波数に対する強さを計算することができます。
最も大きな値を持つ位置を調べるとそこがピークとなり解析した音の周波数になります。
さあ、これを計算するためにMainModuleに次のように書き込んでいきます。
Function EstimatePitch(ByVal FilePath As String, ByVal StartPos As Long, EndPos As Long)
Dim CWAV As New WAV
'データの読み込み
Dim data() As Integer: data = CWAV.Read_WaveFile(FilePath)
Dim Num As Long: Num = 4096 * 2 'FFTのサイズ
Dim InputData() As Double: ReDim InputData(Num) As Double
Dim N, k As Long
Dim x_real() As Double, x_imag() As Double
ReDim x_real(Num) As Double
ReDim x_imag(Num) As Double
Dim Is_Zero As Long: Is_Zero = 0
'波形
For N = StartPos To EndPos
InputData(N - StartPos) = CDbl(data(N))
Is_Zero = Is_Zero + data(N)
Next N
If Is_Zero = 0 Then
EstimatePitch = 0
Exit Function
End If
Call CWAV.FFT(InputData, x_real, x_imag, Num)
'出力
Dim Out() As Double: ReDim Out((Num / 2) - 1) As Double
For k = 0 To (Num / 2) - 1
Out(k) = Sqr((x_real(k)) ^ 2 + (x_imag(k)) ^ 2)
Next k
'データのFFT
Dim Pitch As Long: Pitch = CWAV.Get_Peek_Hz(Out, Num)
EstimatePitch = Pitch
End Function
Sub test()
On Error Resume Next
Dim CWAV As New WAV
Const FilePath As String = "D:\Song.wav"
Dim Pitch As Long, PitchCode As String
Dim i As Long, Pos As Long
Dim cnt As Long: cnt = 16
For i = 0 To cnt * 10 - 1
'短期間で区切って
Pitch = EstimatePitch(FilePath, Pos, Pos + 44100 * 0.1)
PitchCode = CWAV.Get_PitchCode(Pitch)
Cells(i + 1, 5).Value = (i + 1) * 0.1
If PitchCode = "None" Then
Cells(i + 1, 6).Value = 0
Cells(i + 1, 7).Value = PitchCode
Else
Cells(i + 1, 6).Value = Pitch
Cells(i + 1, 7).Value = PitchCode
End If
Pos = Pos + (44100 * 0.1)
Next
End Sub
解析の準備は整いました。読み込むSong.wavのパスは適宜変更してくださいね。
今回周波数と音階の対応関係の登録については説明しませんでしたが、これはダウンロードしていただいたファイルに既に記載されています。
こんなやつです。
では実際にプログラムを走らせてみましょう。testプロシージャを実行します。するとうまく動くと次のようになるはずです。
※一部のみ抜粋
E列が時間で、F列が周波数、G列が推定される音階ですね。カエルの歌は「ドレミファミレド ミファソラソファミ ドドドド ドレミファミレド」なので正しく推定できているようです。こちらをグラフにしてみましょう。
このようになります。いい感じです！ちゃんと音階が正しく読めているようです。
しかしながら今回の内容だけでは和音を含む音や、一般的な曲では解析がうまくできません。興味のある方は色々調べてみると面白いかと思います。
以上で解析は終了です。普通はExcelでフーリエ変換をする機会もないと思いますが、原理の理解とその実装を考える上では機能が少ないが取り扱いがしやすいVBAも有用だと思います。
最後まで読んでくださりありがとうございました。また別の記事でお会いしましょう！
【生涯無期限サポート付き】 Excel/VBAセミナー残席情報
- 【おすすめ】オンライン版利用（動画視聴期限：無期限/会場開催にも参加できます） たった1日で即戦力になるExcelマスター講座
- 【残席9】2023年9月30日（土）13:00-18:00 東京会場 たった1日で即戦力になるExcelマスター講座
- 【残席10】2023年10月7日（土）13:00-18:00 東京会場 たった1日で即戦力になるExcelマスター講座
- 【残席5】2023年10月9日（月）13:00-17:00 福岡会場 ※要PC持参 たった1日で即戦力になるExcelマスター講座
- 【残席10】2023年10月14日（土）13:00-18:00 東京会場 たった1日で即戦力になるExcelマスター講座
- 【残席10】2023年10月21日（土）13:00-18:00 東京会場 たった1日で即戦力になるExcelマスター講座
- 【残席10】2023年11月25日（土）10:00-17:00 東京会場 Excel VBAマスター講座
-
『たった1日で即戦力になるExcelの教科書』執筆陣が自ら直接指導。
実務直結・一日集中・受講後無期限サポート付きのExcelセミナー
