{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# MFCC theory and implementation\n## Theory\nMel Frequency Cepstral Coefficents (MFCCs) is a way of extracting features from an audio. The MFCC uses the MEL scale to divide the frequency band to sub-bands and then extracts the Cepstral Coefficents using Discrete Cosine Transform (DCT). MEL scale is based on the way humans distinguish between frequencies which makes it very convenient to process sounds.\n\nLets first understand how humans perceive sounds.\n\n## Human voice sound perception\nAdult humans fundamental voice frequency range is between 85Hz to 255Hz (85Hz to 180Hz for male and 165Hz to 255Hz for female). On top of the fundamental frequency there are harmonics of fundamental frequencys. Harmonics are whole multiplications of the fundamental frequency. If for instance the fundamental frequency is 100Hz then its second harmonic will be 200Hz, third harmonic is 300Hz and so on.\n\nYou can see an example in the image below [1] which shows frequency vs. time of several pronounced words and color represents frequency power at that point (yellow strongest and black weakest):\n![Speech Spectrogram](https://static.scientificamerican.com/sciam/assets/media/sound/speechsep-audio/speechsep-2-spect.png \"Speech Spectrogram\")\n\nNotice the first horizontal yellow line on the bottom of each segment. That is the fundamental frequency and its the strongest. Above that there are harmonics with the same frequncy distance from each other.\n\nHumans can hear roughly between 20Hz to 20KHz. The perception of sound is non-linear [2] and you can better distinguish between low frequency sounds than high frequency sounds e.g. humans can clearly hear the difference betwee 100Hz and 200Hz but not between 15kHz and 15.1kHz.\n\nYou can try in usin a tone generator: http://www.szynalski.com/tone-generator/\n\n## MEL scale\nA MEL scale is a unit of PITCH proposed by Stevens, Volkmann and Newmann in 1937. The MEL scale is a scale of pitches judged by listeners to be equal in distance one from another [3] [4]. Because of how humans perceive sound the MEL scale is a non-lenear scale and the distances between the pitches increeses with frequency.\n\n# MFCC implementation [5] [6]\n## Sample signal\nFirst, lets load a sample audio and start working with it:","metadata":{"_uuid":"462a26d13ff25773f0e816f55e0355bb4175bbef"}},{"cell_type":"code","source":"import os\nimport numpy as np\nimport scipy\nfrom scipy.io import wavfile\nimport scipy.fftpack as fft\nfrom scipy.signal import get_window\nimport IPython.display as ipd\nimport matplotlib.pyplot as plt\n\n%matplotlib inline","metadata":{"_uuid":"600f9d1247883dc7d47e64e488caf1f29db86baa","collapsed":true,"jupyter":{"outputs_hidden":true}},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"TRAIN_PATH = '../input/audio_train/'\nipd.Audio(TRAIN_PATH + \"a439d172.wav\")","metadata":{"_uuid":"0dfbc8495993481e10e943be56520a6cd56f0218"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_rate, audio = wavfile.read(TRAIN_PATH + \"a439d172.wav\")\nprint(\"Sample rate: {0}Hz\".format(sample_rate))\nprint(\"Audio duration: {0}s\".format(len(audio) / sample_rate))","metadata":{"scrolled":true,"_uuid":"5e68fd22941317512a550e6ad11f1764016e7eed"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Bonus: notice that the sample rate is 44.1kHz which is the standart sampling rate for audio files. This number is not a coincidence. Acording to Nyquist, that sampling rate of 44.1kHz give us band-width of 22.05kHz. The 0Hz to 22.05kHz zone is called the first Nyquist zone, 22.05kHz to 44.1kHz is the second Nyquist zone and so on. If there is a signal in the second Nyquist zone (lets say 22.1kHz) it will alias to the first Nyquist zone (22kHz in our case). We dont want this to happend so we add a low pass filter (LPF). All LPFs has a pass-band, stop-band and transition band. \n\n![LPF Example](https://upload.wikimedia.org/wikipedia/commons/6/60/Butterworth_response.svg \"LPF Example\")\n\nThe sampling frequency of 44.1kHz was chose so that the transition band will be from 20kHz (pass-band) to 22.05kHz (stop-band). It is important to have at least half of the transition band inside the first Nyquist zone because there is stil not enought attenuetion in the filter in this part so there might be aliases from the second Nyquist zone.","metadata":{"_uuid":"28248df2774a2e579d71cf678f071d17bd39b7cb"}},{"cell_type":"code","source":"def normalize_audio(audio):\n    audio = audio / np.max(np.abs(audio))\n    return audio","metadata":{"_uuid":"6c13aa2e69c5cd2b48084d7dce2b00e0f75871e6","collapsed":true,"jupyter":{"outputs_hidden":true}},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"audio = normalize_audio(audio)\nplt.figure(figsize=(15,4))\nplt.plot(np.linspace(0, len(audio) / sample_rate, num=len(audio)), audio)\nplt.grid(True)","metadata":{"_uuid":"9693828e82bafd772fff85dcd71dbf660a2128ac"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There is a lot going on here. It sounds that we should get high and low fundamental frequencies.\n\n## Audio Framing\nBecause audio is a non stationary process, the FFT will produce distortions. To overcome this we can assume that the audio is a stationary process for a short periods of time. Because of that we devide the signal into short frames. Each audio frame will be the same size as the FFT. Also we want the frames to overlap. We do that so that the frames will have some correlation between them and because we loose the information on the edges of each frame after applying a window function.","metadata":{"_uuid":"e94fa032b7e0a5003e7177e70759da3d143fbe83"}},{"cell_type":"code","source":"def frame_audio(audio, FFT_size=2048, hop_size=10, sample_rate=44100):\n    # hop_size in ms\n    \n    audio = np.pad(audio, int(FFT_size / 2), mode='reflect')\n    frame_len = np.round(sample_rate * hop_size / 1000).astype(int)\n    frame_num = int((len(audio) - FFT_size) / frame_len) + 1\n    frames = np.zeros((frame_num,FFT_size))\n    \n    for n in range(frame_num):\n        frames[n] = audio[n*frame_len:n*frame_len+FFT_size]\n    \n    return frames","metadata":{"_uuid":"c40f0f53be2f4923532b4da8fa5bc3130c1ad13e","collapsed":true,"jupyter":{"outputs_hidden":true}},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"hop_size = 15 #ms\nFFT_size = 2048\n\naudio_framed = frame_audio(audio, FFT_size=FFT_size, hop_size=hop_size, sample_rate=sample_rate)\nprint(\"Framed audio shape: {0}\".format(audio_framed.shape))","metadata":{"_uuid":"d85e14f3f0b7f0423cbb58f26b390a9ce92f7920"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we have framed audio matrix with the size of frame number X FFT size.","metadata":{"_uuid":"1c2fe820f0f489ef4fc5ec29baa0670e69444748"}},{"cell_type":"code","source":"print(\"First frame:\")\naudio_framed[1]","metadata":{"_uuid":"860a55d20e132a04e0bc825d6967b1c3834d370a"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Last frame:\")\naudio_framed[-1]","metadata":{"_uuid":"ffacc9ad30af525ec065ec9665baca493af64889"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Convert to frequency domain\nNow we will convert the audio, which is currently in the time domain, to frequency domain. The FFT assumes the audio to be periodic and continues. By framing the signal we assured the audio to be periodic. To make the audio continues, we apply a window function on every frame. If we wont do that, We will get high frequency distortions. To overcome this, we first need to apply a window function to the framed audio and then perforn FFT. The window assures that both ends of the signal will end close to zero.\n\nChoosing the correct window is hard and takes time. For simplicity we will choose the Hanning window. [7]","metadata":{"_uuid":"98585701d4bb3c87fe9cc25ba6ff996dae3c01cb"}},{"cell_type":"code","source":"window = get_window(\"hann\", FFT_size, fftbins=True)\nplt.figure(figsize=(15,4))\nplt.plot(window)\nplt.grid(True)","metadata":{"_uuid":"47534c7e11a9bd3638e447a444675b904bb875c8"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"audio_win = audio_framed * window\n\nind = 69\nplt.figure(figsize=(15,6))\nplt.subplot(2, 1, 1)\nplt.plot(audio_framed[ind])\nplt.title('Original Frame')\nplt.grid(True)\nplt.subplot(2, 1, 2)\nplt.plot(audio_win[ind])\nplt.title('Frame After Windowing')\nplt.grid(True)","metadata":{"_uuid":"ba706661290f58feac29c3fcaeafee3804552077"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In the plot above you can see the that both ends of the frame end on different places on the y axis. You can see that the window brought the edges of each frame closer to zero.\n\nNow lets perform the FFT. After we do the FFT we only take the the positive part of the spectrum (first half +1).","metadata":{"_uuid":"26891a8e507b7b696d9ef7daf5d95ada140eb805"}},{"cell_type":"code","source":"audio_winT = np.transpose(audio_win)\n\naudio_fft = np.empty((int(1 + FFT_size // 2), audio_winT.shape[1]), dtype=np.complex64, order='F')\n\nfor n in range(audio_fft.shape[1]):\n    audio_fft[:, n] = fft.fft(audio_winT[:, n], axis=0)[:audio_fft.shape[0]]\n\naudio_fft = np.transpose(audio_fft)","metadata":{"_uuid":"ff80ef2ff43ecb1b7b99003d1651e122ef0e91be","collapsed":true,"jupyter":{"outputs_hidden":true}},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Calculate signal power","metadata":{"_uuid":"68918ccff758f90428817f112f1c1f790d870cb5"}},{"cell_type":"code","source":"audio_power = np.square(np.abs(audio_fft))\nprint(audio_power.shape)","metadata":{"_uuid":"e4a39b8616c7c9cea613020f1d72a4af2c9f60da"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## MEL-spaced filterbank\nHere we compute the MEL-spaced filterbank and then pass the framed audio through them. That will give us information about the power in each frequency band. The filters can be constructed for any frequency band but for our example we will look on the entire sampled band.\n\nWhat spessial with the MEL-spaced filterbank is the spacing between the filters which grows exponentially with frequency. The filterbank can be made for any frequency band. Here we will compute the filterbank for the entire frequency band.","metadata":{"_uuid":"0ac61886f147ea5de8d06c94ecda5dc24467ae7c"}},{"cell_type":"code","source":"freq_min = 0\nfreq_high = sample_rate / 2\nmel_filter_num = 10\n\nprint(\"Minimum frequency: {0}\".format(freq_min))\nprint(\"Maximum frequency: {0}\".format(freq_high))","metadata":{"_uuid":"16e3fe5d9c2aa4207a4786f83a37633fe94d8660"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Compute filter points\nFirst we construct filter points that determines the start and stop of the filters. To do that we first convert the two filterbank edges to the MEL space. After that we construct a lineary spaced array between the two MEL frequencies. Then we convert the array to the frequency space and finally we normalize the array to the FFT size and choose the associated FFT values.\n\nThis process is shown in the picture below:\n\n<a href=\"https://ibb.co/kCsp4H\"><img src=\"https://preview.ibb.co/hrFnrx/MEL.png\" alt=\"MEL\" border=\"0\" /></a>","metadata":{"_uuid":"204e6a498265d66a3ec402fe7dd4897929420fea"}},{"cell_type":"code","source":"def freq_to_mel(freq):\n    return 2595.0 * np.log10(1.0 + freq / 700.0)\n\ndef met_to_freq(mels):\n    return 700.0 * (10.0**(mels / 2595.0) - 1.0)","metadata":{"_uuid":"923bcf6c33898bc99068a52236a981d22ad2625d","collapsed":true,"jupyter":{"outputs_hidden":true}},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_filter_points(fmin, fmax, mel_filter_num, FFT_size, sample_rate=44100):\n    fmin_mel = freq_to_mel(fmin)\n    fmax_mel = freq_to_mel(fmax)\n    \n    print(\"MEL min: {0}\".format(fmin_mel))\n    print(\"MEL max: {0}\".format(fmax_mel))\n    \n    mels = np.linspace(fmin_mel, fmax_mel, num=mel_filter_num+2)\n    freqs = met_to_freq(mels)\n    \n    return np.floor((FFT_size + 1) / sample_rate * freqs).astype(int), freqs","metadata":{"_uuid":"6c258837e70b0c6a88e32c82cf9d2c31dfed8756","collapsed":true,"jupyter":{"outputs_hidden":true}},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"filter_points, mel_freqs = get_filter_points(freq_min, freq_high, mel_filter_num, FFT_size, sample_rate=44100)\nfilter_points","metadata":{"_uuid":"b5a7fcd45faacf93605e2f438ba881047eb2144e"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Construct the filterbank\nAfter we have the filter points, we construct the filters.","metadata":{"_uuid":"f3585fd874c338aa0f52552e0f48d9137a9ee163"}},{"cell_type":"code","source":"def get_filters(filter_points, FFT_size):\n    filters = np.zeros((len(filter_points)-2,int(FFT_size/2+1)))\n    \n    for n in range(len(filter_points)-2):\n        filters[n, filter_points[n] : filter_points[n + 1]] = np.linspace(0, 1, filter_points[n + 1] - filter_points[n])\n        filters[n, filter_points[n + 1] : filter_points[n + 2]] = np.linspace(1, 0, filter_points[n + 2] - filter_points[n + 1])\n    \n    return filters","metadata":{"_uuid":"904dd238008a4eb82a350b9fed582b751c65654f","collapsed":true,"jupyter":{"outputs_hidden":true}},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"filters = get_filters(filter_points, FFT_size)\n\nplt.figure(figsize=(15,4))\nfor n in range(filters.shape[0]):\n    plt.plot(filters[n])","metadata":{"scrolled":true,"_uuid":"ec545047665b553c980d25c74e09861df9aa94e5"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next we divide the triangular MEL weights by the width of the MEL band (area normalization). If we wont normalize the filters, we will see the noise increase with frequency because of the filter width.","metadata":{"_uuid":"51d1ddf8cc20c534b50d3e7dd9fff74094e04b38"}},{"cell_type":"code","source":"# taken from the librosa library\nenorm = 2.0 / (mel_freqs[2:mel_filter_num+2] - mel_freqs[:mel_filter_num])\nfilters *= enorm[:, np.newaxis]","metadata":{"_uuid":"e84d8c18f9df476b35313d920adf610913967941","collapsed":true,"jupyter":{"outputs_hidden":true}},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15,4))\nfor n in range(filters.shape[0]):\n    plt.plot(filters[n])","metadata":{"_uuid":"1d4af49cae387b6ee2fb6165946e9b72a88a847c"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Filter the signal","metadata":{"_uuid":"b2459afe24f27349434ae429c4f996442613ec18"}},{"cell_type":"code","source":"audio_filtered = np.dot(filters, np.transpose(audio_power))\naudio_log = 10.0 * np.log10(audio_filtered)\naudio_log.shape","metadata":{"_uuid":"977f1e8ff6aee33644aa419a4130191ac79e3005"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we have a matrix represemting the audio power in all 10 filters in different time frames.\n## Generate the Cepstral Coefficents\nThe final step in generating the MFCC is to use the Discrete Cosine Transform (DCT). We will use the DCT-III. This type of DCT will extract high frequency and low frequency changes in the the signal [8] [9].","metadata":{"_uuid":"e5927662e8f9248c1bf8c0668af5d306ea31d620"}},{"cell_type":"code","source":"def dct(dct_filter_num, filter_len):\n    basis = np.empty((dct_filter_num,filter_len))\n    basis[0, :] = 1.0 / np.sqrt(filter_len)\n    \n    samples = np.arange(1, 2 * filter_len, 2) * np.pi / (2.0 * filter_len)\n\n    for i in range(1, dct_filter_num):\n        basis[i, :] = np.cos(i * samples) * np.sqrt(2.0 / filter_len)\n        \n    return basis","metadata":{"_uuid":"e4815ea74fad45820bbcd2a3417296235277d915","collapsed":true,"jupyter":{"outputs_hidden":true}},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dct_filter_num = 40\n\ndct_filters = dct(dct_filter_num, mel_filter_num)\n\ncepstral_coefficents = np.dot(dct_filters, audio_log)\ncepstral_coefficents.shape","metadata":{"_uuid":"b922616cfc75cf76a8e42570e09397726ce8c244"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Reviewing Cepstral coefficents","metadata":{"_uuid":"68fd684fba3c5fa004997108d3ed24bc8d1e4e45"}},{"cell_type":"code","source":"cepstral_coefficents[:, 0]","metadata":{"_uuid":"9e8363263103be23b810783de1b6bc33de9f361b"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15,5))\nplt.plot(np.linspace(0, len(audio) / sample_rate, num=len(audio)), audio)\nplt.imshow(cepstral_coefficents, aspect='auto', origin='lower');","metadata":{"_uuid":"6a302215afb0e64c65b6f2b87f91c6941a690cb9"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Conclution\nThe MFCC is a good algorithm to extract high frequency and low frequency information. There are many parts that can be changed depending on the application. One thing that we might want to check is how the power is scattered across frequency band and choosing the filter amount based on that.\n\n# Bibliography\n[1] https://en.wikipedia.org/wiki/Voice_frequency\n\n[2] https://en.wikipedia.org/wiki/Hearing_range\n\n[3] https://www.sfu.ca/sonic-studio/handbook/Mel.html\n\n[4] https://en.wikipedia.org/wiki/Mel_scale\n\n[5] http://practicalcryptography.com/miscellaneous/machine-learning/guide-mel-frequency-cepstral-coefficients-mfccs/\n\n[6] https://librosa.github.io/librosa/\n\n[7] https://en.wikipedia.org/wiki/Window_function#Hamming_window\n\n[8] https://www.youtube.com/watch?v=Q2aEzeMDHMA\n\n[9] https://en.wikipedia.org/wiki/Discrete_cosine_transform","metadata":{"_uuid":"02441af327f4a298f6e4ef6b35d6d6694f0cbc32"}}]}