获取声音特征MFCC VB.NET

Imports System.Numerics
Imports MathNet.Numerics.IntegralTransforms

Public Class MFCC


    Public H As Double(,)

    Private MFCCNum As Integer
    Private FrameSize As Integer '帧长512


    Public Sub New(Optional framesize As Integer = 512, Optional MFCCNum As Integer = 26)
        '注意设置最小频率 freMin  0  ,300 
        Me.MFCCNum = MFCCNum
        Me.FrameSize = framesize

        H = New Double(MFCCNum, Me.FrameSize / 2) {}

        '计算mel系数
        Dim filter_points(40 + 1) As Integer '40个滤波器，需要41点

        Const sampleRate As Integer = 16000  '采样频率 16000
        Const filterNum As Integer = 40  '滤波器数量 取40个

        Dim freMax As Double = sampleRate / 2   '实际最大频率 
        Dim freMin As Double = 0    '实际最小频率 
        Dim melFremax As Double = 1125 * Math.Log(1 + freMax / 700)     '将实际频率转换成梅尔频率 
        Dim melFremin As Double = 1125 * Math.Log(1 + freMin / 700)


        Dim k As Double = (melFremax - melFremin) / (filterNum + 1)


        Dim m As Double() = New Double(filterNum + 1) {}
        Dim r As Double() = New Double(filterNum + 1) {}


        For i As Integer = 0 To filterNum + 1
            m(i) = melFremin + k * i
            r(i) = 700 * (Math.Exp(m(i) / 1125) - 1)
            '将梅尔频率转换成实际频率 
            filter_points(i) = Math.Floor((Me.FrameSize + 1) * r(i) / sampleRate)
        Next


        '生成mel滤波器
        For i As Integer = 0 To MFCCNum
            For j As Integer = 0 To Me.FrameSize / 2 - 1
                If j < filter_points(i) Then
                    H(i, j) = 0
                End If
                If (filter_points(i) <= j) And (j <= filter_points(i + 1)) Then
                    H(i, j) = (CDbl(j - filter_points(i)) / (filter_points(i + 1) - filter_points(i)))
                End If
                If (filter_points(i + 1) <= j) And (j <= filter_points(i + 2)) Then
                    H(i, j) = (CDbl(filter_points(i + 2) - j) / (filter_points(i + 2) - filter_points(i + 1)))
                End If
                If j > filter_points(i + 2) Then
                    H(i, j) = 0
                End If
            Next
        Next



    End Sub

    '汉明窗
    Public Sub Hamming_window(WaveData() As Single)
        Dim len As Integer = WaveData.Length
        Dim omega As Single = 2.0 * Math.PI / len
        For j As Integer = 0 To len - 1
            WaveData(j) = (0.54 - 0.46 * Math.Cos(omega * (j))) * WaveData(j)
        Next
    End Sub

    '傅里叶计算
    Public Function FFT(WaveData() As Single) As Complex()
        Dim FFT_Complex(WaveData.Length - 1) As Complex
        For i = 0 To WaveData.Length - 1
            FFT_Complex(i) = WaveData(i)
        Next
        MathNet.Numerics.IntegralTransforms.Fourier.Forward(FFT_Complex, FourierOptions.Matlab)
        Return FFT_Complex


    End Function


    Public Function MFCC(fft() As Complex) As Single()

        '取LOG
        Dim S As Single() = New Single(MFCCNum - 1) {}
        For i As Integer = 0 To MFCCNum - 1
            For j As Integer = 0 To Me.FrameSize / 2 - 1
                S(i) = S(i) + Math.Pow(fft(j).Magnitude, 2) * H(i, j)
            Next
            If S(i) <> 0 Then
                S(i) = Math.Log(S(i), Math.E)
            End If
        Next


        'DCT运算
        Dim mfcc_mass(MFCCNum - 1) As Double
        For l As Integer = 0 To MFCCNum - 1
            For i As Integer = 0 To MFCCNum - 1
                mfcc_mass(l) += S(i) * Math.Cos(Math.PI * l * ((i * 0.5) / 20))




            Next
        Next

        Return S


    End Function

End Class

使用方法：

 Public Function WavTMfcc(data() As Int16) As Single()
        '分帧，每一帧进行mfcc计算 帧长512 帧移256 
        Dim len = data.Length
        Dim FrmSize = 512
        Dim FrmNum = len / 256 - 1

        Dim mfccs(FrmNum * 26 - 1) As Single  '7*26 =182
        Dim Frame(FrmSize - 1) As Single
        For i As Integer = 0 To FrmNum - 1
            Array.Copy(data, i * 256, Frame, 0, FrmSize)
            mfcc.Hamming_window(Frame)
            Dim fft As Complex() = mfcc.FFT(Frame)
            Dim rs As Single() = mfcc.MFCC(fft)
            Array.Copy(rs, 0, mfccs, i * 26, 26)

        Next

        Return mfccs
    End Function