mel frequency cepstral coefficients mfcc calculation mfcc is an algorithm widely used in audio and speech processing to represent the shortterm power spectrum of a sound signal in a more compact and discriminative way it is particularly popular in speech and audio processing tasks such as speech recognition and speaker identification how mel frequency cepstral coefficients are calculated 1 preprocessing load an audio signal and normalize it to ensure that the values fall within a specific range e g between 1 and 1 frame the audio signal into overlapping fixedlength segments typically using a technique like windowing to reduce spectral leakage 2 fourier transform apply a fast fourier transform fft to each audio frame to convert it from the time domain to the frequency domain this results in a representation of the audio frame as a sequence of frequency components 3 power spectrum calculate the power spectrum by taking the squared magnitude of each frequency component obtained from the fft this step measures the energy distribution across different frequency bands 4 mel filterbank apply a set of triangular filterbanks spaced in the mel frequency scale to the power spectrum these filters mimic the human auditory system s frequency response each filterbank sums the power spectrum values within its band 5 logarithmic compression take the logarithm typically base 10 of the filterbank values to compress the dynamic range this step mimics the logarithmic response of the human ear to sound intensity 6 discrete cosine transform dct apply the discrete cosine transform to the log filterbank energies to obtain the mfcc coefficients this transformation helps decorrelate the filterbank energies and captures the most important features of the audio signal 7 feature extraction select a subset of the dct coefficients to form the feature vector often the first few coefficients e g 1213 are used for most applications references melfrequency cepstral coefficients mfccs https en wikipedia orgwikimelfrequencycepstrum speech and language processing by daniel jurafsky james h martin https web stanford edujurafskyslp3 mel frequency cepstral coefficient mfcc tutorial http practicalcryptography commiscellaneousmachinelearning guidemelfrequencycepstralcoefficientsmfccs amir lavasani calculate mel frequency cepstral coefficients mfccs from an audio signal args audio the input audio signal samplerate the sample rate of the audio signal in hz fttsize the size of the fft window default is 1024 hoplength the hop length for frame creation default is 20ms melfilternum the number of mel filters default is 10 dctfilternum the number of dct filters default is 40 returns a matrix of mfccs for the input audio raises valueerror if the input audio is empty example samplerate 44100 sample rate of 44 1 khz duration 2 0 duration of 1 second t np linspace0 duration intsamplerate duration endpointfalse audio 0 5 np sin2 np pi 440 0 t generate a 440 hz sine wave mfccs mfccaudio samplerate mfccs shape 40 101 normalize audio frame audio into convert to frequency domain for simplicity we will choose the hanning window normalize an audio signal by scaling it to have values between 1 and 1 args audio the input audio signal returns the normalized audio signal examples audio np array1 2 3 4 5 normalizedaudio normalizeaudio np maxnormalizedaudio 1 0 np minnormalizedaudio 0 2 divide the entire audio signal by the maximum absolute value split an audio signal into overlapping frames args audio the input audio signal samplerate the sample rate of the audio signal hoplength the length of the hopping default is 20ms fttsize the size of the fft window default is 1024 returns an array of overlapping frames examples audio np array1 2 3 4 5 6 7 8 9 101000 samplerate 8000 frames audioframesaudio samplerate hoplength10 fttsize512 frames shape 126 512 pad the audio signal to handle edge cases calculate the number of frames initialize an array to store the frames split the audio signal into frames calculate the fast fourier transform fft of windowed audio data args audiowindowed the windowed audio signal fttsize the size of the fft default is 1024 returns the fft of the audio data examples audiowindowed np array1 0 2 0 3 0 4 0 5 0 6 0 audiofft calculatefftaudiowindowed fttsize4 np allcloseaudiofft0 np array6 00 j 1 50 8660254j 1 50 8660254j true transpose the audio data to have time in rows and channels in columns initialize an array to store the fft results compute fft for each channel transpose the fft results back to the original shape calculate the power of the audio signal from its fft args audiofft the fft of the audio signal returns the power of the audio signal examples audiofft np array12j 23j 34j 45j power calculatesignalpoweraudiofft np allclosepower np array5 13 25 41 true calculate the power by squaring the absolute values of the fft coefficients convert a frequency in hertz to the mel scale args freq the frequency in hertz returns the frequency in mel scale examples roundfreqtomel1000 2 999 99 use the formula to convert frequency to the mel scale convert a frequency in the mel scale to hertz args mels the frequency in mel scale returns the frequency in hertz examples roundmeltofreq999 99 2 1000 01 use the formula to convert mel scale to frequency create a melspaced filter bank for audio processing args samplerate the sample rate of the audio melfilternum the number of mel filters default is 10 fttsize the size of the fft default is 1024 returns melspaced filter bank examples roundmelspacedfilterbank8000 10 102401 10 0 0004603981 calculate filter points and mel frequencies normalize filters taken from the librosa library generate filters for audio processing args filterpoints a list of filter points fttsize the size of the fft returns a matrix of filters examples getfiltersnp array0 20 51 95 161 256 dtypeint 512 shape 4 257 linearly increase values from 0 to 1 linearly decrease values from 1 to 0 calculate the filter points and frequencies for mel frequency filters args samplerate the sample rate of the audio freqmin the minimum frequency in hertz freqhigh the maximum frequency in hertz melfilternum the number of mel filters default is 10 fttsize the size of the fft default is 1024 returns filter points and corresponding frequencies examples filterpoints getfilterpoints8000 0 4000 melfilternum4 fttsize512 filterpoints0 array 0 20 51 95 161 256 filterpoints1 array 0 324 46707094 799 33254207 1494 30973963 2511 42581671 4000 convert minimum and maximum frequencies to mel scale generate equally spaced mel frequencies convert mel frequencies back to hertz calculate filter points as integer values compute the discrete cosine transform dct basis matrix args dctfilternum the number of dct filters to generate filternum the number of the fbank filters returns the dct basis matrix examples rounddiscretecosinetransform3 500 5 0 44721 example function to calculate mel frequency cepstral coefficients mfccs from an audio file args wavfilepath the path to the wav audio file returns np ndarray the computed mfccs for the audio load the audio from the wav file calculate mfccs calculate mel frequency cepstral coefficients mfccs from an audio signal args audio the input audio signal sample_rate the sample rate of the audio signal in hz ftt_size the size of the fft window default is 1024 hop_length the hop length for frame creation default is 20ms mel_filter_num the number of mel filters default is 10 dct_filter_num the number of dct filters default is 40 returns a matrix of mfccs for the input audio raises valueerror if the input audio is empty example sample_rate 44100 sample rate of 44 1 khz duration 2 0 duration of 1 second t np linspace 0 duration int sample_rate duration endpoint false audio 0 5 np sin 2 np pi 440 0 t generate a 440 hz sine wave mfccs mfcc audio sample_rate mfccs shape 40 101 normalize audio frame audio into convert to frequency domain for simplicity we will choose the hanning window normalize an audio signal by scaling it to have values between 1 and 1 args audio the input audio signal returns the normalized audio signal examples audio np array 1 2 3 4 5 normalized_audio normalize audio np max normalized_audio 1 0 np min normalized_audio 0 2 divide the entire audio signal by the maximum absolute value split an audio signal into overlapping frames args audio the input audio signal sample_rate the sample rate of the audio signal hop_length the length of the hopping default is 20ms ftt_size the size of the fft window default is 1024 returns an array of overlapping frames examples audio np array 1 2 3 4 5 6 7 8 9 10 1000 sample_rate 8000 frames audio_frames audio sample_rate hop_length 10 ftt_size 512 frames shape 126 512 pad the audio signal to handle edge cases calculate the number of frames initialize an array to store the frames split the audio signal into frames calculate the fast fourier transform fft of windowed audio data args audio_windowed the windowed audio signal ftt_size the size of the fft default is 1024 returns the fft of the audio data examples audio_windowed np array 1 0 2 0 3 0 4 0 5 0 6 0 audio_fft calculate_fft audio_windowed ftt_size 4 np allclose audio_fft 0 np array 6 0 0 j 1 5 0 8660254j 1 5 0 8660254j true transpose the audio data to have time in rows and channels in columns initialize an array to store the fft results compute fft for each channel transpose the fft results back to the original shape calculate the power of the audio signal from its fft args audio_fft the fft of the audio signal returns the power of the audio signal examples audio_fft np array 1 2j 2 3j 3 4j 4 5j power calculate_signal_power audio_fft np allclose power np array 5 13 25 41 true calculate the power by squaring the absolute values of the fft coefficients convert a frequency in hertz to the mel scale args freq the frequency in hertz returns the frequency in mel scale examples round freq_to_mel 1000 2 999 99 use the formula to convert frequency to the mel scale convert a frequency in the mel scale to hertz args mels the frequency in mel scale returns the frequency in hertz examples round mel_to_freq 999 99 2 1000 01 use the formula to convert mel scale to frequency create a mel spaced filter bank for audio processing args sample_rate the sample rate of the audio mel_filter_num the number of mel filters default is 10 ftt_size the size of the fft default is 1024 returns mel spaced filter bank examples round mel_spaced_filterbank 8000 10 1024 0 1 10 0 0004603981 calculate filter points and mel frequencies normalize filters taken from the librosa library generate filters for audio processing args filter_points a list of filter points ftt_size the size of the fft returns a matrix of filters examples get_filters np array 0 20 51 95 161 256 dtype int 512 shape 4 257 linearly increase values from 0 to 1 linearly decrease values from 1 to 0 calculate the filter points and frequencies for mel frequency filters args sample_rate the sample rate of the audio freq_min the minimum frequency in hertz freq_high the maximum frequency in hertz mel_filter_num the number of mel filters default is 10 ftt_size the size of the fft default is 1024 returns filter points and corresponding frequencies examples filter_points get_filter_points 8000 0 4000 mel_filter_num 4 ftt_size 512 filter_points 0 array 0 20 51 95 161 256 filter_points 1 array 0 324 46707094 799 33254207 1494 30973963 2511 42581671 4000 convert minimum and maximum frequencies to mel scale generate equally spaced mel frequencies convert mel frequencies back to hertz calculate filter points as integer values compute the discrete cosine transform dct basis matrix args dct_filter_num the number of dct filters to generate filter_num the number of the fbank filters returns the dct basis matrix examples round discrete_cosine_transform 3 5 0 0 5 0 44721 example function to calculate mel frequency cepstral coefficients mfccs from an audio file args wav_file_path the path to the wav audio file returns np ndarray the computed mfccs for the audio load the audio from the wav file calculate mfccs
import logging

import numpy as np
import scipy.fftpack as fft
from scipy.signal import get_window

logging.basicConfig(filename=f"{__file__}.log", level=logging.INFO)


def mfcc(
    audio: np.ndarray,
    sample_rate: int,
    ftt_size: int = 1024,
    hop_length: int = 20,
    mel_filter_num: int = 10,
    dct_filter_num: int = 40,
) -> np.ndarray:
    
    logging.info(f"Sample rate: {sample_rate}Hz")
    logging.info(f"Audio duration: {len(audio) / sample_rate}s")
    logging.info(f"Audio min: {np.min(audio)}")
    logging.info(f"Audio max: {np.max(audio)}")

    
    audio_normalized = normalize(audio)

    logging.info(f"Normalized audio min: {np.min(audio_normalized)}")
    logging.info(f"Normalized audio max: {np.max(audio_normalized)}")

    
    audio_framed = audio_frames(
        audio_normalized, sample_rate, ftt_size=ftt_size, hop_length=hop_length
    )

    logging.info(f"Framed audio shape: {audio_framed.shape}")
    logging.info(f"First frame: {audio_framed[0]}")

    
    
    window = get_window("hann", ftt_size, fftbins=True)
    audio_windowed = audio_framed * window

    logging.info(f"Windowed audio shape: {audio_windowed.shape}")
    logging.info(f"First frame: {audio_windowed[0]}")

    audio_fft = calculate_fft(audio_windowed, ftt_size)
    logging.info(f"fft audio shape: {audio_fft.shape}")
    logging.info(f"First frame: {audio_fft[0]}")

    audio_power = calculate_signal_power(audio_fft)
    logging.info(f"power audio shape: {audio_power.shape}")
    logging.info(f"First frame: {audio_power[0]}")

    filters = mel_spaced_filterbank(sample_rate, mel_filter_num, ftt_size)
    logging.info(f"filters shape: {filters.shape}")

    audio_filtered = np.dot(filters, np.transpose(audio_power))
    audio_log = 10.0 * np.log10(audio_filtered)
    logging.info(f"audio_log shape: {audio_log.shape}")

    dct_filters = discrete_cosine_transform(dct_filter_num, mel_filter_num)
    cepstral_coefficents = np.dot(dct_filters, audio_log)

    logging.info(f"cepstral_coefficents shape: {cepstral_coefficents.shape}")
    return cepstral_coefficents


def normalize(audio: np.ndarray) -> np.ndarray:
    
    
    return audio / np.max(np.abs(audio))


def audio_frames(
    audio: np.ndarray,
    sample_rate: int,
    hop_length: int = 20,
    ftt_size: int = 1024,
) -> np.ndarray:
    

    hop_size = np.round(sample_rate * hop_length / 1000).astype(int)

    
    audio = np.pad(audio, int(ftt_size / 2), mode="reflect")

    
    frame_count = int((len(audio) - ftt_size) / hop_size) + 1

    
    frames = np.zeros((frame_count, ftt_size))

    
    for n in range(frame_count):
        frames[n] = audio[n * hop_size : n * hop_size + ftt_size]

    return frames


def calculate_fft(audio_windowed: np.ndarray, ftt_size: int = 1024) -> np.ndarray:
    
    
    audio_transposed = np.transpose(audio_windowed)

    
    audio_fft = np.empty(
        (int(1 + ftt_size // 2), audio_transposed.shape[1]),
        dtype=np.complex64,
        order="F",
    )

    
    for n in range(audio_fft.shape[1]):
        audio_fft[:, n] = fft.fft(audio_transposed[:, n], axis=0)[: audio_fft.shape[0]]

    
    return np.transpose(audio_fft)


def calculate_signal_power(audio_fft: np.ndarray) -> np.ndarray:
    
    
    return np.square(np.abs(audio_fft))


def freq_to_mel(freq: float) -> float:
    
    
    return 2595.0 * np.log10(1.0 + freq / 700.0)


def mel_to_freq(mels: float) -> float:
    
    
    return 700.0 * (10.0 ** (mels / 2595.0) - 1.0)


def mel_spaced_filterbank(
    sample_rate: int, mel_filter_num: int = 10, ftt_size: int = 1024
) -> np.ndarray:
    
    freq_min = 0
    freq_high = sample_rate // 2

    logging.info(f"Minimum frequency: {freq_min}")
    logging.info(f"Maximum frequency: {freq_high}")

    
    filter_points, mel_freqs = get_filter_points(
        sample_rate,
        freq_min,
        freq_high,
        mel_filter_num,
        ftt_size,
    )

    filters = get_filters(filter_points, ftt_size)

    
    
    enorm = 2.0 / (mel_freqs[2 : mel_filter_num + 2] - mel_freqs[:mel_filter_num])
    return filters * enorm[:, np.newaxis]


def get_filters(filter_points: np.ndarray, ftt_size: int) -> np.ndarray:
    
    num_filters = len(filter_points) - 2
    filters = np.zeros((num_filters, int(ftt_size / 2) + 1))

    for n in range(num_filters):
        start = filter_points[n]
        mid = filter_points[n + 1]
        end = filter_points[n + 2]

        
        filters[n, start:mid] = np.linspace(0, 1, mid - start)

        
        filters[n, mid:end] = np.linspace(1, 0, end - mid)

    return filters


def get_filter_points(
    sample_rate: int,
    freq_min: int,
    freq_high: int,
    mel_filter_num: int = 10,
    ftt_size: int = 1024,
) -> tuple[np.ndarray, np.ndarray]:
    
    
    fmin_mel = freq_to_mel(freq_min)
    fmax_mel = freq_to_mel(freq_high)

    logging.info(f"MEL min: {fmin_mel}")
    logging.info(f"MEL max: {fmax_mel}")

    
    mels = np.linspace(fmin_mel, fmax_mel, num=mel_filter_num + 2)

    
    freqs = mel_to_freq(mels)

    
    filter_points = np.floor((ftt_size + 1) / sample_rate * freqs).astype(int)

    return filter_points, freqs


def discrete_cosine_transform(dct_filter_num: int, filter_num: int) -> np.ndarray:
    
    basis = np.empty((dct_filter_num, filter_num))
    basis[0, :] = 1.0 / np.sqrt(filter_num)

    samples = np.arange(1, 2 * filter_num, 2) * np.pi / (2.0 * filter_num)

    for i in range(1, dct_filter_num):
        basis[i, :] = np.cos(i * samples) * np.sqrt(2.0 / filter_num)

    return basis


def example(wav_file_path: str = "./path-to-file/sample.wav") -> np.ndarray:
    
    from scipy.io import wavfile

    
    sample_rate, audio = wavfile.read(wav_file_path)

    
    return mfcc(audio, sample_rate)


if __name__ == "__main__":
    import doctest

    doctest.testmod()