{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":8900,"databundleVersionId":862232,"sourceType":"competition"}],"dockerImageVersionId":30746,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# MFCC theory\n\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","metadata":{}},{"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":{"execution":{"iopub.status.busy":"2024-07-18T05:13:16.000643Z","iopub.execute_input":"2024-07-18T05:13:16.001280Z","iopub.status.idle":"2024-07-18T05:13:16.008622Z","shell.execute_reply.started":"2024-07-18T05:13:16.001246Z","shell.execute_reply":"2024-07-18T05:13:16.007517Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"TRAIN_PATH = '/kaggle/input/freesound-audio-tagging/audio_train/'\nsample_rate, audio = wavfile.read(TRAIN_PATH + \"a439d172.wav\")\nipd.Audio(TRAIN_PATH + \"a439d172.wav\" , rate = sample_rate)","metadata":{"execution":{"iopub.status.busy":"2024-07-18T05:15:15.824431Z","iopub.execute_input":"2024-07-18T05:15:15.824936Z","iopub.status.idle":"2024-07-18T05:15:15.877558Z","shell.execute_reply.started":"2024-07-18T05:15:15.824890Z","shell.execute_reply":"2024-07-18T05:15:15.875617Z"},"trusted":true},"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":{"execution":{"iopub.status.busy":"2024-07-18T05:15:45.150053Z","iopub.execute_input":"2024-07-18T05:15:45.150427Z","iopub.status.idle":"2024-07-18T05:15:45.169606Z","shell.execute_reply.started":"2024-07-18T05:15:45.150397Z","shell.execute_reply":"2024-07-18T05:15:45.168564Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def normalize_audio(audio):\n    audio = audio / np.max(np.abs(audio))\n    return audio","metadata":{"execution":{"iopub.status.busy":"2024-07-18T05:16:13.559474Z","iopub.execute_input":"2024-07-18T05:16:13.560368Z","iopub.status.idle":"2024-07-18T05:16:13.565134Z","shell.execute_reply.started":"2024-07-18T05:16:13.560335Z","shell.execute_reply":"2024-07-18T05:16:13.564042Z"},"trusted":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":{"execution":{"iopub.status.busy":"2024-07-18T05:16:26.490321Z","iopub.execute_input":"2024-07-18T05:16:26.491403Z","iopub.status.idle":"2024-07-18T05:16:26.992638Z","shell.execute_reply.started":"2024-07-18T05:16:26.491368Z","shell.execute_reply":"2024-07-18T05:16:26.991412Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Audio Framing","metadata":{}},{"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":{"execution":{"iopub.status.busy":"2024-07-18T05:16:43.997985Z","iopub.execute_input":"2024-07-18T05:16:43.998391Z","iopub.status.idle":"2024-07-18T05:16:44.006320Z","shell.execute_reply.started":"2024-07-18T05:16:43.998357Z","shell.execute_reply":"2024-07-18T05:16:44.005149Z"},"trusted":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":{"execution":{"iopub.status.busy":"2024-07-18T05:16:51.742094Z","iopub.execute_input":"2024-07-18T05:16:51.742871Z","iopub.status.idle":"2024-07-18T05:16:51.752044Z","shell.execute_reply.started":"2024-07-18T05:16:51.742839Z","shell.execute_reply":"2024-07-18T05:16:51.750945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"First frame:\")\naudio_framed[1]","metadata":{"execution":{"iopub.status.busy":"2024-07-18T05:17:01.160206Z","iopub.execute_input":"2024-07-18T05:17:01.160638Z","iopub.status.idle":"2024-07-18T05:17:01.169391Z","shell.execute_reply.started":"2024-07-18T05:17:01.160604Z","shell.execute_reply":"2024-07-18T05:17:01.168138Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Last frame:\")\naudio_framed[-1]","metadata":{"execution":{"iopub.status.busy":"2024-07-18T05:17:07.119764Z","iopub.execute_input":"2024-07-18T05:17:07.120156Z","iopub.status.idle":"2024-07-18T05:17:07.128168Z","shell.execute_reply.started":"2024-07-18T05:17:07.120128Z","shell.execute_reply":"2024-07-18T05:17:07.127195Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Convert to frequency domain","metadata":{}},{"cell_type":"code","source":"window = get_window(\"hann\", FFT_size, fftbins=True)\nplt.figure(figsize=(15,4))\nplt.plot(window)\nplt.grid(True)","metadata":{"execution":{"iopub.status.busy":"2024-07-18T05:17:25.800681Z","iopub.execute_input":"2024-07-18T05:17:25.801146Z","iopub.status.idle":"2024-07-18T05:17:26.073876Z","shell.execute_reply.started":"2024-07-18T05:17:25.801110Z","shell.execute_reply":"2024-07-18T05:17:26.072757Z"},"trusted":true},"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":{"execution":{"iopub.status.busy":"2024-07-18T05:17:31.699621Z","iopub.execute_input":"2024-07-18T05:17:31.699984Z","iopub.status.idle":"2024-07-18T05:17:32.215254Z","shell.execute_reply.started":"2024-07-18T05:17:31.699956Z","shell.execute_reply":"2024-07-18T05:17:32.214162Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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":{"execution":{"iopub.status.busy":"2024-07-18T05:17:37.666381Z","iopub.execute_input":"2024-07-18T05:17:37.667447Z","iopub.status.idle":"2024-07-18T05:17:37.694551Z","shell.execute_reply.started":"2024-07-18T05:17:37.667409Z","shell.execute_reply":"2024-07-18T05:17:37.693566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Calculate signal power","metadata":{}},{"cell_type":"code","source":"audio_power = np.square(np.abs(audio_fft))\nprint(audio_power.shape)","metadata":{"execution":{"iopub.status.busy":"2024-07-18T05:17:55.149994Z","iopub.execute_input":"2024-07-18T05:17:55.150407Z","iopub.status.idle":"2024-07-18T05:17:55.157035Z","shell.execute_reply.started":"2024-07-18T05:17:55.150375Z","shell.execute_reply":"2024-07-18T05:17:55.155840Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## MEL-spaced filterbank","metadata":{}},{"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":{"execution":{"iopub.status.busy":"2024-07-18T05:18:12.109708Z","iopub.execute_input":"2024-07-18T05:18:12.110536Z","iopub.status.idle":"2024-07-18T05:18:12.116750Z","shell.execute_reply.started":"2024-07-18T05:18:12.110487Z","shell.execute_reply":"2024-07-18T05:18:12.115389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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":{"execution":{"iopub.status.busy":"2024-07-18T05:18:30.100137Z","iopub.execute_input":"2024-07-18T05:18:30.100583Z","iopub.status.idle":"2024-07-18T05:18:30.106507Z","shell.execute_reply.started":"2024-07-18T05:18:30.100549Z","shell.execute_reply":"2024-07-18T05:18:30.105358Z"},"trusted":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":{"execution":{"iopub.status.busy":"2024-07-18T05:18:34.872661Z","iopub.execute_input":"2024-07-18T05:18:34.873115Z","iopub.status.idle":"2024-07-18T05:18:34.880430Z","shell.execute_reply.started":"2024-07-18T05:18:34.873084Z","shell.execute_reply":"2024-07-18T05:18:34.879239Z"},"trusted":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":{"execution":{"iopub.status.busy":"2024-07-18T05:18:42.103384Z","iopub.execute_input":"2024-07-18T05:18:42.104569Z","iopub.status.idle":"2024-07-18T05:18:42.114195Z","shell.execute_reply.started":"2024-07-18T05:18:42.104530Z","shell.execute_reply":"2024-07-18T05:18:42.112981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Construct the filterbank","metadata":{}},{"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":{"execution":{"iopub.status.busy":"2024-07-18T05:19:05.659598Z","iopub.execute_input":"2024-07-18T05:19:05.660289Z","iopub.status.idle":"2024-07-18T05:19:05.667892Z","shell.execute_reply.started":"2024-07-18T05:19:05.660255Z","shell.execute_reply":"2024-07-18T05:19:05.666628Z"},"trusted":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":{"execution":{"iopub.status.busy":"2024-07-18T05:19:10.019839Z","iopub.execute_input":"2024-07-18T05:19:10.020242Z","iopub.status.idle":"2024-07-18T05:19:10.275404Z","shell.execute_reply.started":"2024-07-18T05:19:10.020201Z","shell.execute_reply":"2024-07-18T05:19:10.274332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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":{"execution":{"iopub.status.busy":"2024-07-18T05:19:15.262593Z","iopub.execute_input":"2024-07-18T05:19:15.263512Z","iopub.status.idle":"2024-07-18T05:19:15.269197Z","shell.execute_reply.started":"2024-07-18T05:19:15.263451Z","shell.execute_reply":"2024-07-18T05:19:15.268044Z"},"trusted":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":{"execution":{"iopub.status.busy":"2024-07-18T05:19:21.379809Z","iopub.execute_input":"2024-07-18T05:19:21.380302Z","iopub.status.idle":"2024-07-18T05:19:21.687033Z","shell.execute_reply.started":"2024-07-18T05:19:21.380262Z","shell.execute_reply":"2024-07-18T05:19:21.685683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Filter the signal","metadata":{}},{"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":{"execution":{"iopub.status.busy":"2024-07-18T05:19:38.991660Z","iopub.execute_input":"2024-07-18T05:19:38.992562Z","iopub.status.idle":"2024-07-18T05:19:39.009847Z","shell.execute_reply.started":"2024-07-18T05:19:38.992525Z","shell.execute_reply":"2024-07-18T05:19:39.008385Z"},"trusted":true},"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","metadata":{}},{"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":{"execution":{"iopub.status.busy":"2024-07-18T05:19:59.000790Z","iopub.execute_input":"2024-07-18T05:19:59.001605Z","iopub.status.idle":"2024-07-18T05:19:59.008708Z","shell.execute_reply.started":"2024-07-18T05:19:59.001571Z","shell.execute_reply":"2024-07-18T05:19:59.007388Z"},"trusted":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":{"execution":{"iopub.status.busy":"2024-07-18T05:20:03.829407Z","iopub.execute_input":"2024-07-18T05:20:03.830417Z","iopub.status.idle":"2024-07-18T05:20:03.840098Z","shell.execute_reply.started":"2024-07-18T05:20:03.830383Z","shell.execute_reply":"2024-07-18T05:20:03.838244Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Reviewing Cepstral coefficents","metadata":{}},{"cell_type":"code","source":"cepstral_coefficents[:, 0]","metadata":{"execution":{"iopub.status.busy":"2024-07-18T05:20:18.839486Z","iopub.execute_input":"2024-07-18T05:20:18.839866Z","iopub.status.idle":"2024-07-18T05:20:18.847454Z","shell.execute_reply.started":"2024-07-18T05:20:18.839835Z","shell.execute_reply":"2024-07-18T05:20:18.846268Z"},"trusted":true},"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":{"execution":{"iopub.status.busy":"2024-07-18T05:20:24.460015Z","iopub.execute_input":"2024-07-18T05:20:24.460676Z","iopub.status.idle":"2024-07-18T05:20:25.026072Z","shell.execute_reply.started":"2024-07-18T05:20:24.460642Z","shell.execute_reply":"2024-07-18T05:20:25.024871Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}