Excelの効率化がわかる記事一覧


Excelで簡易的にWAVファイルから周波数分析をしてみよう

こんにちは。すごい改善のAitoです。

今日は少し変わったネタで書いていきます。それも「ExcelWAVファイルから周波数分析を行う」です!一体何を言っているんだ?という方も、これやったことあるという方も様々な反応が出ると思います。理系で学んでいると専攻によってはどこかで音声処理を学ぶ機会があります。ここでは最も身近な分析ソフトである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セミナー残席情報

  • 【おすすめ】オンライン版利用(動画視聴期限:無期限/会場開催にも参加できます)
  • 【残席5】2024年2月23日(祝)13:00-18:00 東京会場 たった1日で即戦力になるExcelマスター講座
  • 【残席10】2024年3月10日(日)10:00-17:00 東京会場 面倒なExcel作業を自動化するマクロVBAマスター講座
  • 【残席10】2024年3月16日(土)13:00-18:00 東京会場 たった1日で即戦力になるExcelマスター講座
  • 【残席10】2024年3月23日(土)13:00-18:00 東京会場 たった1日で即戦力になるExcelマスター講座
  • 【残席10】2024年3月30日(土)13:00-18:00 東京会場 たった1日で即戦力になるExcelマスター講座
  • 『たった1日で即戦力になるExcelの教科書』執筆陣が自ら直接指導。
    実務直結・一日集中・受講後無期限サポート付きのExcelセミナー

法人向けExcelオンライン研修

利用人数無制限・個別サポート付きのExcelオンライン研修です。

すごい改善のExcelセミナー

1日完結・実務直結・無期限サポートつきのExcelセミナー開催中

累計40万部突破のExcel教科書

Excel関数Tシャツ等のご購入はこちら