{"cells":[{"metadata":{"_uuid":"cf808b6b464476e44c8dda4beada9884cf91adca","_cell_guid":"cec0e376-8e36-4751-8b53-2f298fda3a47"},"cell_type":"markdown","source":"## Introduction \n\nHello! This is my very first Kernel. It is meant to give a grasp of a problem of speech representation. I'd also like to take a look on a features specific to this dataset. \n\nContent:<br>\n* [1. Visualization of the recordings - input features](#visualization)\n   * [1.1. Wave and spectrogram](#waveandspectrogram)\n   * [1.2. MFCC](#mfcc)\n   * [1.3. Sprectrogram in 3d](#3d)\n   * [1.4. Silence removal](#resampl)\n   * [1.5. Resampling - dimensionality reductions](#silenceremoval)\n   * [1.6. Features extraction steps](#featuresextractionsteps)\n* [2. Dataset investigation](#investigations)\n   * [2.1. Number of files](#numberoffiles)\n   * [2.2. Mean spectrograms and fft](#meanspectrogramsandfft)\n   * [2.3. Deeper into recordings](#deeper)\n   * [2.4. Length of recordings](#len)\n   * [2.5. Note on Gaussian Mixtures modeling](#gmms)\n   * [2.6. Frequency components across the words](#components)\n   * [2.7. Anomaly detection](#anomaly)\n* [3. Where to look for the inspiration](#wheretostart)\n\nAll we need is here:"},{"metadata":{"collapsed":true,"_uuid":"d9596d80cb6445d4214dda15e40d777cadbd4669","_cell_guid":"8fd82027-7be0-4a4e-a921-b8acacaaf077","trusted":false},"cell_type":"code","source":"import os\nfrom os.path import isdir, join\nfrom pathlib import Path\nimport pandas as pd\n\n# Math\nimport numpy as np\nfrom scipy.fftpack import fft\nfrom scipy import signal\nfrom scipy.io import wavfile\nimport librosa\n\nfrom sklearn.decomposition import PCA\n\n# Visualization\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport IPython.display as ipd\nimport librosa.display\n\nimport plotly.offline as py\npy.init_notebook_mode(connected=True)\nimport plotly.graph_objs as go\nimport plotly.tools as tls\nimport pandas as pd\n\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"95fabaca63ab1a486bcc1f6824b26919ef325ff4","_cell_guid":"7f050711-6810-4aac-a306-da82fddb5579"},"cell_type":"markdown","source":"\n# 1. Visualization \n<a id=\"visualization\"></a> \n\nThere are two theories of a human hearing - place ( https://en.wikipedia.org/wiki/Place_theory_(hearing) (frequency-based) and temporal (https://en.wikipedia.org/wiki/Temporal_theory_(hearing) )\nIn speech recognition, I see two main tendencies - to input [spectrogram](https://en.wikipedia.org/wiki/Spectrogram) (frequencies), and more sophisticated features MFCC - Mel-Frequency Cepstral Coefficients, PLP. You rarely work with raw, temporal data.\n\nLet's visualize some recordings!\n\n## 1.1. Wave and spectrogram:\n<a id=\"waveandspectrogram\"></a> \n\nChoose and read some file:"},{"metadata":{"collapsed":true,"_uuid":"76266716e7df45a83073fb2964218c85b36d31cb","_cell_guid":"02126a6d-dd84-4f0a-88eb-ed9ff46a9bdf","trusted":false},"cell_type":"code","source":"train_audio_path = '../input/train/audio/'\nfilename = '/yes/0a7c2a8d_nohash_0.wav'\nsample_rate, samples = wavfile.read(str(train_audio_path) + filename)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3bc26d76ea9f627c4d476ff8e9523f37d0668bbf","_cell_guid":"a7715152-3866-48dd-8bbb-31a72e9aa9bf"},"cell_type":"markdown","source":"Define a function that calculates spectrogram.\n\nNote, that we are taking logarithm of spectrogram values. It will make our plot much more clear, moreover, it is strictly connected to the way people hear.\nWe need to assure that there are no 0 values as input to logarithm."},{"metadata":{"collapsed":true,"_uuid":"a3569f66d5bbbdcf338eaa121328a507f3a7b431","_cell_guid":"e464fe63-138e-4c66-a1f7-3ad3a81daa38","trusted":false},"cell_type":"code","source":"def log_specgram(audio, sample_rate, window_size=20,\n                 step_size=10, eps=1e-10):\n    nperseg = int(round(window_size * sample_rate / 1e3))\n    noverlap = int(round(step_size * sample_rate / 1e3))\n    freqs, times, spec = signal.spectrogram(audio,\n                                    fs=sample_rate,\n                                    window='hann',\n                                    nperseg=nperseg,\n                                    noverlap=noverlap,\n                                    detrend=False)\n    return freqs, times, np.log(spec.T.astype(np.float32) + eps)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4fd53946fd96b09765a267231ea5a66b313c2d4e","_cell_guid":"625dcb59-00ec-4b3f-97d5-f8adc12ac61a"},"cell_type":"markdown","source":"Frequencies are in range (0, 8000) according to [Nyquist theorem](https://en.wikipedia.org/wiki/Nyquist_rate).\n\nLet's plot it:"},{"metadata":{"collapsed":true,"_uuid":"ec1d065704c51f7f2b5b49f00da64809257815ab","_cell_guid":"4f77267a-1720-439b-9ef9-90e60f4446e1","trusted":false},"cell_type":"code","source":"freqs, times, spectrogram = log_specgram(samples, sample_rate)\n\nfig = plt.figure(figsize=(14, 8))\nax1 = fig.add_subplot(211)\nax1.set_title('Raw wave of ' + filename)\nax1.set_ylabel('Amplitude')\nax1.plot(np.linspace(0, sample_rate/len(samples), sample_rate), samples)\n\nax2 = fig.add_subplot(212)\nax2.imshow(spectrogram.T, aspect='auto', origin='lower', \n           extent=[times.min(), times.max(), freqs.min(), freqs.max()])\nax2.set_yticks(freqs[::16])\nax2.set_xticks(times[::16])\nax2.set_title('Spectrogram of ' + filename)\nax2.set_ylabel('Freqs in Hz')\nax2.set_xlabel('Seconds')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f36fd74c9ad998d71b3a8838347b1bdbe8c82a7","_cell_guid":"013846a9-a929-45d9-97f5-98c59c6b2f23"},"cell_type":"markdown","source":"If we use spectrogram as an input features for NN, we have to remember to normalize features. (We need to normalize over all the dataset, here's example just for one, which doesn't give good *mean* and *std*!)"},{"metadata":{"collapsed":true,"_uuid":"b3d09cf8bd1e91f54774f84dc952508d8a8a4eb8","_cell_guid":"9572b5e1-0b0f-42aa-934e-a1e313c21f46","trusted":false},"cell_type":"code","source":"mean = np.mean(spectrogram, axis=0)\nstd = np.std(spectrogram, axis=0)\nspectrogram = (spectrogram - mean) / std","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ee5c8c77122ece10687bacd92cd7ec01ab77af75","_cell_guid":"e3acd63a-4e3d-48fc-ba01-c28fedcd4496"},"cell_type":"markdown","source":"There is an interesting fact to point out. We have ~160 features for each frame, frequencies are between 0 and 8000. It means, that one feature corresponds to 50 Hz. However, [frequency resolution of the ear is 3.6 Hz within the octave of 1000 – 2000 Hz](https://en.wikipedia.org/wiki/Psychoacoustics) It means, that people are far more precise and can hear much smaller details than those represented by spectrograms like above."},{"metadata":{"_uuid":"4eb99845d61397b9acb2488d34e2bafa7aad4cca","_cell_guid":"53904969-d453-4f0e-8e9b-6d932190bed1"},"cell_type":"markdown","source":"## 1.2. MFCC\n<a id=\"mfcc\"></a> \n\nIf you want to get to know some details about *MFCC* take a look at this great tutorial. [MFCC explained](http://practicalcryptography.com/miscellaneous/machine-learning/guide-mel-frequency-cepstral-coefficients-mfccs/) You can see, that it is well prepared to imitate human hearing properties.\n\nYou can calculate *Mel power spectrogram* and *MFCC* using for example *librosa* python package.\n"},{"metadata":{"collapsed":true,"_uuid":"4d996e6499140446685a3796418faa15a5f9d425","_cell_guid":"a6cb80ed-0e64-43b5-87ae-d33f3f844276","trusted":false},"cell_type":"code","source":"# From this tutorial\n# https://github.com/librosa/librosa/blob/master/examples/LibROSA%20demo.ipynb\nS = librosa.feature.melspectrogram(samples, sr=sample_rate, n_mels=128)\n\n# Convert to log scale (dB). We'll use the peak power (max) as reference.\nlog_S = librosa.power_to_db(S, ref=np.max)\n\nplt.figure(figsize=(12, 4))\nlibrosa.display.specshow(log_S, sr=sample_rate, x_axis='time', y_axis='mel')\nplt.title('Mel power spectrogram ')\nplt.colorbar(format='%+02.0f dB')\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_uuid":"e1a21ec3fdbb30b360479d8886e3e496a5511ba4","_cell_guid":"38c436c0-9db0-48f8-a3d2-8ad660447bea","trusted":false},"cell_type":"code","source":"mfcc = librosa.feature.mfcc(S=log_S, n_mfcc=13)\n\n# Let's pad on the first and second deltas while we're at it\ndelta2_mfcc = librosa.feature.delta(mfcc, order=2)\n\nplt.figure(figsize=(12, 4))\nlibrosa.display.specshow(delta2_mfcc)\nplt.ylabel('MFCC coeffs')\nplt.xlabel('Time')\nplt.title('MFCC')\nplt.colorbar()\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ea023dd3edc2aea6a82b20eb4c53aac7f818390e","_cell_guid":"d1db710e-20d5-40ee-ac57-2e1e5de19c0e"},"cell_type":"markdown","source":"In classical, but still state-of-the-art systems, *MFCC*  or similar features are taken as the input to the system instead of spectrograms.\n\nHowever, in end-to-end (often neural-network based) systems, the most common input features are probably raw spectrograms, or mel power spectrograms. For example *MFCC* decorrelates features, but NNs deal with correlated features well. Also, if you'll understand mel filters, you may consider their usage sensible.a\n\nIt is your decision which to choose!"},{"metadata":{"_uuid":"62f174047f35264a71dad240503d319a452271f9","_cell_guid":"156af766-21f8-4f93-ba1a-750157a11171"},"cell_type":"markdown","source":"## 1.3. Spectrogram in 3d\n<a id=\"3d\"></a> \n\nBy the way, times change, and the tools change. Have you ever seen spectrogram in 3d?"},{"metadata":{"collapsed":true,"_uuid":"cfb8b4827cecd2e01a1416176035488edd454fe4","_cell_guid":"51a5bc43-5216-4f13-8bc5-f84e005a01df","trusted":false},"cell_type":"code","source":"data = [go.Surface(z=spectrogram.T)]\nlayout = go.Layout(\n    title='Specgtrogram of \"yes\" in 3d',\n    scene = dict(\n    yaxis = dict(title='Frequencies', range=freqs),\n    xaxis = dict(title='Time', range=times),\n    zaxis = dict(title='Log amplitude'),\n    ),\n)\nfig = go.Figure(data=data, layout=layout)\npy.iplot(fig)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"790e79ac09b6cea0b509965a86760bdcbe2671a7","_cell_guid":"7b8a7587-f41a-4021-9694-637f29e2c3f8"},"cell_type":"markdown","source":"(Don't know how to set axis ranges to proper values yet. I'd also like it to be streched like a classic spectrogram above..)"},{"metadata":{"_uuid":"769e6738c4dae9923b9c0b0a99981bce8b443030","_cell_guid":"a2ad2019-f402-4226-9bad-65fb400aa8b1"},"cell_type":"markdown","source":"## 1.4. Silence removal\n<a id=\"silenceremoval\"></a> \n\nLet's listen to that file"},{"metadata":{"collapsed":true,"_uuid":"ab0145dc0c8efdc08b4153b136c2b78634f6ed07","scrolled":false,"_cell_guid":"f49b916e-53a2-4dbe-bd03-3d8d93bf25a6","trusted":false},"cell_type":"code","source":"ipd.Audio(samples, rate=sample_rate)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9745fb19ce26c85c312a20e7fa19d98e672ceb64","_cell_guid":"4c23e8b3-0c8f-4eda-8f35-7486bdecfd9d"},"cell_type":"markdown","source":"I consider that some *VAD* (Voice Activity Detection) will be really useful here. Although the words are short, there is a lot of silence in them. A decent *VAD* can reduce training size a lot, accelerating training speed significantly.\nLet's cut a bit of the file from the beginning and from the end. and listen to it again (based on a plot above, we take from 4000 to 13000):"},{"metadata":{"collapsed":true,"_uuid":"539d123d84ac0181b820cca82c4098ab0ca54116","_cell_guid":"2c85e04d-cbd8-4702-bd50-7340497e800d","trusted":false},"cell_type":"code","source":"samples_cut = samples[4000:13000]\nipd.Audio(samples_cut, rate=sample_rate)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e9ceecbabecc11f789b3de382ee4c909186e6d22","_cell_guid":"45898236-4528-4e21-86dd-55abcf4f639f"},"cell_type":"markdown","source":"We can agree that the entire word can be heard. It is impossible to cut all the files manually and do this basing on the simple plot. But you can use for example *webrtcvad* package to have a good *VAD*.\n\nLet's plot it again, together with guessed alignment of* 'y' 'e' 's'* graphems"},{"metadata":{"collapsed":true,"_uuid":"6831fa9311397dc8bca4192f657767d36c5c1a38","_cell_guid":"038fe488-4f25-42bd-af11-108f5ecbb1e7","trusted":false},"cell_type":"code","source":"freqs, times, spectrogram_cut = log_specgram(samples_cut, sample_rate)\n\nfig = plt.figure(figsize=(14, 8))\nax1 = fig.add_subplot(211)\nax1.set_title('Raw wave of ' + filename)\nax1.set_ylabel('Amplitude')\nax1.plot(samples_cut)\n\nax2 = fig.add_subplot(212)\nax2.set_title('Spectrogram of ' + filename)\nax2.set_ylabel('Frequencies * 0.1')\nax2.set_xlabel('Samples')\nax2.imshow(spectrogram_cut.T, aspect='auto', origin='lower', \n           extent=[times.min(), times.max(), freqs.min(), freqs.max()])\nax2.set_yticks(freqs[::16])\nax2.set_xticks(times[::16])\nax2.text(0.06, 1000, 'Y', fontsize=18)\nax2.text(0.17, 1000, 'E', fontsize=18)\nax2.text(0.36, 1000, 'S', fontsize=18)\n\nxcoords = [0.025, 0.11, 0.23, 0.49]\nfor xc in xcoords:\n    ax1.axvline(x=xc*16000, c='r')\n    ax2.axvline(x=xc, c='r')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e8f5fa497bbd2b3f5e7dbb9fa20d59d9773309a1","_cell_guid":"f081f185-336a-429d-ba71-c0d2337c35ae"},"cell_type":"markdown","source":"## 1.5. Resampling - dimensionality reduction\n<a id=\"resampl\"></a> \n\nAnother way to reduce the dimensionality of our data is to resample recordings.\n\nYou can hear that the recording don't sound very natural, because they are sampled with 16k frequency, and we usually hear much more. However, [the most speech related frequencies are presented in smaller band](https://en.wikipedia.org/wiki/Voice_frequency). That's why you can still understand another person talking to the telephone, where GSM signal is sampled to 8000 Hz.\n\nSummarizing, we could resample our dataset to 8k. We will discard some information that shouldn't be important, and we'll reduce size of the data.\n\nWe have to remember that it can be risky, because this is a competition, and sometimes very small difference in performance wins, so we don't want to lost anything. On the other hand, first experiments can be done much faster with smaller training size.\n\nWe'll need to calculate FFT (Fast Fourier Transform). Definition:\n"},{"metadata":{"collapsed":true,"_uuid":"213de6a783443d118c3509acc26f9f4bd0319d85","_cell_guid":"86dedd69-e084-403e-9c02-5370018acf1c","trusted":false},"cell_type":"code","source":"def custom_fft(y, fs):\n    T = 1.0 / fs\n    N = y.shape[0]\n    yf = fft(y)\n    xf = np.linspace(0.0, 1.0/(2.0*T), N//2)\n    vals = 2.0/N * np.abs(yf[0:N//2])  # FFT is simmetrical, so we take just the first half\n    # FFT is also complex, to we take just the real part (abs)\n    return xf, vals","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"665e57b4652493e6d3b61ba2b7e70967170e7900","_cell_guid":"0fc3b446-d19e-4cd2-b1d6-3cf58ff332bf"},"cell_type":"markdown","source":"Let's read some recording, resample it, and listen. We can also compare FFT, Notice, that there is almost no information above 4000 Hz in original signal."},{"metadata":{"collapsed":true,"_uuid":"b8fdb36dc4fce089ea5a3c3dcc27f65625232e34","_cell_guid":"919e85ca-7769-4214-a1d7-5eaa74a32b19","trusted":false},"cell_type":"code","source":"filename = '/happy/0b09edd3_nohash_0.wav'\nnew_sample_rate = 8000\n\nsample_rate, samples = wavfile.read(str(train_audio_path) + filename)\nresampled = signal.resample(samples, int(new_sample_rate/sample_rate * samples.shape[0]))","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_uuid":"afa8138a2ae7888ade44713fb5f8451f9c9e7f02","_cell_guid":"13f397f1-cd5d-4f0f-846a-0edd9f58bcff","trusted":false},"cell_type":"code","source":"ipd.Audio(samples, rate=sample_rate)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_uuid":"3f600c9414ab5cef205c814ba16a356d4121790b","_cell_guid":"5ab11b21-9528-47fa-8ff0-244b1d0c94b3","trusted":false},"cell_type":"code","source":"ipd.Audio(resampled, rate=new_sample_rate)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"96380594085d818693b959307d371e95f727f03b","_cell_guid":"37da8174-e6aa-463d-bef7-c8b20c6ca513"},"cell_type":"markdown","source":"Almost no difference!"},{"metadata":{"collapsed":true,"_uuid":"4448038dfa22ec582cde229346cb1ba309c76b9f","_cell_guid":"baed6102-3c75-4f16-85d7-723d8a084b9a","trusted":false},"cell_type":"code","source":"xf, vals = custom_fft(samples, sample_rate)\nplt.figure(figsize=(12, 4))\nplt.title('FFT of recording sampled with ' + str(sample_rate) + ' Hz')\nplt.plot(xf, vals)\nplt.xlabel('Frequency')\nplt.grid()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_uuid":"88953237ea59d13e9647813bef06a911f06f0e61","_cell_guid":"3cc1a49a-4cd4-49ed-83c8-f2437062f8be","trusted":false},"cell_type":"code","source":"xf, vals = custom_fft(resampled, new_sample_rate)\nplt.figure(figsize=(12, 4))\nplt.title('FFT of recording sampled with ' + str(new_sample_rate) + ' Hz')\nplt.plot(xf, vals)\nplt.xlabel('Frequency')\nplt.grid()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"152c1b14d7a7b57d7ab4fb0bd52e38564406cb92","_cell_guid":"592ffc6a-edda-4b08-9419-d3462599da5c"},"cell_type":"markdown","source":"This is how we reduced dataset size twice!"},{"metadata":{"_uuid":"57fe8c6a25753e2eb46285bc8d725d20182c1421","_cell_guid":"f98fe35d-2d56-4153-b054-0882bd2e58ce"},"cell_type":"markdown","source":"## 1.6. Features extraction steps\n<a id=\"featuresextractionsteps\"></a> \n\nI would propose the feature extraction algorithm like that:\n1. Resampling\n2. *VAD*\n3. Maybe padding with 0 to make signals be equal length\n4. Log spectrogram (or *MFCC*, or *PLP*)\n5. Features normalization with *mean* and *std*\n6. Stacking of a given number of frames to get temporal information\n\nIt's a pity it can't be done in notebook. It has not much sense to write things from zero, and everything is ready to take, but in packages, that can not be imported in Kernels."},{"metadata":{"_uuid":"caf345ca07983f1e1d4f8a05f6f74859554289db","_cell_guid":"3d36bac6-eb6f-4a53-b148-805493e39052"},"cell_type":"markdown","source":"\n# 2. Dataset investigation\n<a id=\"investigations\"></a> \n\nSome usuall investgation of dataset.\n\n## 2.1. Number of records\n<a id=\"numberoffiles\"></a> \n"},{"metadata":{"collapsed":true,"_uuid":"59826a3eb0f60439d5beee06781193bc67cc53f7","_cell_guid":"3c24fbdd-e50e-47a1-8c44-1894bec7f043","trusted":false},"cell_type":"code","source":"dirs = [f for f in os.listdir(train_audio_path) if isdir(join(train_audio_path, f))]\ndirs.sort()\nprint('Number of labels: ' + str(len(dirs)))","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_uuid":"ea30edeaf3d8020bf55ee2a57af230bded9732e2","_cell_guid":"a6b82ced-df8c-4c7a-8d4c-ed32bf9f60f6","trusted":false},"cell_type":"code","source":"# Calculate\nnumber_of_recordings = []\nfor direct in dirs:\n    waves = [f for f in os.listdir(join(train_audio_path, direct)) if f.endswith('.wav')]\n    number_of_recordings.append(len(waves))\n\n# Plot\ndata = [go.Histogram(x=dirs, y=number_of_recordings)]\ntrace = go.Bar(\n    x=dirs,\n    y=number_of_recordings,\n    marker=dict(color = number_of_recordings, colorscale='Viridius', showscale=True\n    ),\n)\nlayout = go.Layout(\n    title='Number of recordings in given label',\n    xaxis = dict(title='Words'),\n    yaxis = dict(title='Number of recordings')\n)\npy.iplot(go.Figure(data=[trace], layout=layout))\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"928d2933df5e0a37c7dc40f2ec50b9d10423d533","_cell_guid":"c4b1d377-d84e-4601-97fa-2c989023c400"},"cell_type":"markdown","source":"Dataset is balanced except of background_noise, but that's the different thing."},{"metadata":{"_uuid":"0847bb69f23ce4c8153f869222c03731fd62a22e","_cell_guid":"b9b6b43c-96fa-489d-8a49-b77dabd41705"},"cell_type":"markdown","source":"## 2.2. Deeper into recordings\n<a id=\"deeper\"></a> "},{"metadata":{"_uuid":"d4b8b90afc03493f93babb7b9d401ebd0caa1c18","_cell_guid":"766c2e18-2aff-43c1-9672-bd51d4348867"},"cell_type":"markdown","source":"There's a very important fact. Recordings come from very different sources. As far as I can tell, some of them can come from mobile GSM channel.\n\nNevertheless,** it is extremely important to split the dataset in a way that one speaker doesn't occur in both train and test sets.**\nJust take a look and listen to this two examlpes:"},{"metadata":{"collapsed":true,"_uuid":"de451cfb747953d60d9b6982c046b417fdc1ab9b","_cell_guid":"d29ea00a-54eb-474a-8770-fbd269bbd21e","trusted":false},"cell_type":"code","source":"filenames = ['on/004ae714_nohash_0.wav', 'on/0137b3f4_nohash_0.wav']\nfor filename in filenames:\n    sample_rate, samples = wavfile.read(str(train_audio_path) + filename)\n    xf, vals = custom_fft(samples, sample_rate)\n    plt.figure(figsize=(12, 4))\n    plt.title('FFT of speaker ' + filename[4:11])\n    plt.plot(xf, vals)\n    plt.xlabel('Frequency')\n    plt.grid()\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ae779083e33d29a22bcf972542d9911a7b9d64de","_cell_guid":"508c79b0-6beb-43cb-9092-9e0c5719e12e"},"cell_type":"markdown","source":"Even better to listen:"},{"metadata":{"collapsed":true,"_uuid":"4e2d4b6d2c6e6806a2b0b0d4554e49b080003a62","_cell_guid":"63e5851f-b462-4f9c-8e17-a7e8e301d616","trusted":false},"cell_type":"code","source":"print('Speaker ' + filenames[0][4:11])\nipd.Audio(join(train_audio_path, filenames[0]))","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_uuid":"c335e31977a721831149579da9105cc2b664b40d","_cell_guid":"fd2eb518-df30-44a0-89b3-9c4611beb50d","trusted":false},"cell_type":"code","source":"print('Speaker ' + filenames[1][4:11])\nipd.Audio(join(train_audio_path, filenames[1]))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"63555f55e906409e008d9c0a988a03b26cbb8983","_cell_guid":"7c2f7df0-d062-46d6-9507-5ab548db14bb"},"cell_type":"markdown","source":"There are also recordings with some weird silence (some compression?):\n"},{"metadata":{"collapsed":true,"_uuid":"ba2e76ebe90eb9c4fed6617d1948e0c71d88c54e","_cell_guid":"bc6075ee-fc43-4d04-b400-64896ac8450d","trusted":false},"cell_type":"code","source":"filename = '/yes/01bb6a2a_nohash_1.wav'\nsample_rate, samples = wavfile.read(str(train_audio_path) + filename)\nfreqs, times, spectrogram = log_specgram(samples, sample_rate)\n\nplt.figure(figsize=(10, 7))\nplt.title('Spectrogram of ' + filename)\nplt.ylabel('Freqs')\nplt.xlabel('Time')\nplt.imshow(spectrogram.T, aspect='auto', origin='lower', \n           extent=[times.min(), times.max(), freqs.min(), freqs.max()])\nplt.yticks(freqs[::16])\nplt.xticks(times[::16])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"cbf303188e33f325691aa9447f18586788bf110a","_cell_guid":"fb57dc50-0510-4dcf-ba07-7999daf7349e"},"cell_type":"markdown","source":"It means, that we have to prevent overfitting to the very specific acoustical environments.\n"},{"metadata":{"_uuid":"85c9f5b8e2dac9bf2aa2edf0ebebc0ae53ff6533","_cell_guid":"ae626131-069f-4243-b7fc-c0014b11e2d8"},"cell_type":"markdown","source":"## 2.3. Recordings length\n<a id=\"len\"></a> \n\nFind if all the files have 1 second duration:"},{"metadata":{"collapsed":true,"_uuid":"16a2e2c908235a99f64024abab272c65d3d99c65","_cell_guid":"23be7e5e-e4b4-40a0-b9a3-4bc850571a28","trusted":false},"cell_type":"code","source":"num_of_shorter = 0\nfor direct in dirs:\n    waves = [f for f in os.listdir(join(train_audio_path, direct)) if f.endswith('.wav')]\n    for wav in waves:\n        sample_rate, samples = wavfile.read(train_audio_path + direct + '/' + wav)\n        if samples.shape[0] < sample_rate:\n            num_of_shorter += 1\nprint('Number of recordings shorter than 1 second: ' + str(num_of_shorter))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c035e161e9c9622fa96f9589ffbfe826e01c5658","_cell_guid":"57b26071-e603-4ee7-b1d5-5bd2bd46e438"},"cell_type":"markdown","source":"That's suprising, and there is a lot of them. We can pad them with zeros."},{"metadata":{"_uuid":"96655097f9173d34d529a0626446194473cebf69","_cell_guid":"771ff180-73dd-4be8-aa56-ccecb3586416"},"cell_type":"markdown","source":"## 2.4. Mean spectrograms and FFT\n<a id=\"meanspectrogramsandfft\"></a> "},{"metadata":{"_uuid":"a3f64232afa284102e8ffdcdbe0db509f4a78a7e","_cell_guid":"a621a03b-4812-4d7c-93dd-6b0bf7f10572"},"cell_type":"markdown","source":"Let's plot mean FFT for every word"},{"metadata":{"collapsed":true,"_uuid":"e8fc86e267bb48d3ae35f8e2d85db070f28889c9","_cell_guid":"dc065096-6888-4a38-9c66-0729bfa858f6","trusted":false},"cell_type":"code","source":"to_keep = 'yes no up down left right on off stop go'.split()\ndirs = [d for d in dirs if d in to_keep]\n\nprint(dirs)\n\nfor direct in dirs:\n    vals_all = []\n    spec_all = []\n\n    waves = [f for f in os.listdir(join(train_audio_path, direct)) if f.endswith('.wav')]\n    for wav in waves:\n        sample_rate, samples = wavfile.read(train_audio_path + direct + '/' + wav)\n        if samples.shape[0] != 16000:\n            continue\n        xf, vals = custom_fft(samples, 16000)\n        vals_all.append(vals)\n        freqs, times, spec = log_specgram(samples, 16000)\n        spec_all.append(spec)\n\n    plt.figure(figsize=(14, 4))\n    plt.subplot(121)\n    plt.title('Mean fft of ' + direct)\n    plt.plot(np.mean(np.array(vals_all), axis=0))\n    plt.grid()\n    plt.subplot(122)\n    plt.title('Mean specgram of ' + direct)\n    plt.imshow(np.mean(np.array(spec_all), axis=0).T, aspect='auto', origin='lower', \n               extent=[times.min(), times.max(), freqs.min(), freqs.max()])\n    plt.yticks(freqs[::16])\n    plt.xticks(times[::16])\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"931aceb2fb6ac23defc699b3d423b510171b1626","_cell_guid":"1089a473-7fe5-40a0-9e43-8f3cbf1901c3"},"cell_type":"markdown","source":"## 2.5. Gaussian Mixtures modeling\n<a id=\"gmms\"></a> \n\nWe can see that mean FFT looks different for every word. We could model each FFT with a mixture of Gaussian distributions. Some of them however, look almost identical on FFT, like *stop* and *up*... But wait, they are still distinguishable when we look at spectrograms! High frequencies are earlier than low at the beginning of *stop* (probably *s*).\n\nThat's why temporal component is also necessary. There is a [Kaldi](http://kaldi-asr.org/) library, that can model words (or smaller parts of words) with GMMs and model temporal dependencies with [Hidden Markov Models](https://github.com/danijel3/ASRDemos/blob/master/notebooks/HMM_FST.ipynb).\n\nWe could use simple GMMs for words to check what can we model and how hard it is to distinguish the words. We can use [Scikit-learn](http://scikit-learn.org/) for that, however it is not straightforward and lasts very long here, so I abandon this idea for now."},{"metadata":{"_uuid":"0c89774ecfd33f29c10cba58f1fb1c12647b0928","_cell_guid":"341fd72e-750e-4d78-a7b7-1d347bcad4e8"},"cell_type":"markdown","source":"## 2.6. Frequency components across the words\n<a id=\"components\"></a> \n"},{"metadata":{"collapsed":true,"_uuid":"0c3e4ecb1eb87bd8c86547df68db35ef9a84a2ce","_cell_guid":"7a883de7-1fa7-428d-bd92-0fb0d3f75aed","trusted":false},"cell_type":"code","source":"def violinplot_frequency(dirs, freq_ind):\n    \"\"\" Plot violinplots for given words (waves in dirs) and frequency freq_ind\n    from all frequencies freqs.\"\"\"\n\n    spec_all = []  # Contain spectrograms\n    ind = 0\n    for direct in dirs:\n        spec_all.append([])\n\n        waves = [f for f in os.listdir(join(train_audio_path, direct)) if\n                 f.endswith('.wav')]\n        for wav in waves[:100]:\n            sample_rate, samples = wavfile.read(\n                train_audio_path + direct + '/' + wav)\n            freqs, times, spec = log_specgram(samples, sample_rate)\n            spec_all[ind].extend(spec[:, freq_ind])\n        ind += 1\n\n    # Different lengths = different num of frames. Make number equal\n    minimum = min([len(spec) for spec in spec_all])\n    spec_all = np.array([spec[:minimum] for spec in spec_all])\n\n    plt.figure(figsize=(13,7))\n    plt.title('Frequency ' + str(freqs[freq_ind]) + ' Hz')\n    plt.ylabel('Amount of frequency in a word')\n    plt.xlabel('Words')\n    sns.violinplot(data=pd.DataFrame(spec_all.T, columns=dirs))\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_uuid":"8932034e13e12d7ff50bc6663c03153ae1052378","_cell_guid":"08aae780-53a3-4047-b548-46b8893aaed9","trusted":false},"cell_type":"code","source":"violinplot_frequency(dirs, 20)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_uuid":"30e1bf01aa6431fd883530beabb4da5d9c8518dc","_cell_guid":"ea7ea51e-3a38-46db-891d-4c3aef8fd810","trusted":false},"cell_type":"code","source":"violinplot_frequency(dirs, 50)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_uuid":"fb2a917be56bc19c256ad52e257122a7f1b5dd8b","_cell_guid":"a84e722a-3aaa-42dd-868b-6eb629abf40d","trusted":false},"cell_type":"code","source":"violinplot_frequency(dirs, 120)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f10445f9398dcd57b404591c69075a23ab3d5115","_cell_guid":"7e06a571-5ca6-4a50-8bf1-c037d9399df3"},"cell_type":"markdown","source":"## 2.7. Anomaly detection\n<a id=\"anomaly\"></a> \n\nWe should check if there are any recordings that somehow stand out from the rest. We can lower the dimensionality of the dataset and interactively check for any anomaly.\nWe'll use PCA for dimensionality reduction:"},{"metadata":{"collapsed":true,"_uuid":"5a52c1fbea1b68c78502ac9954bd7ecae19149d6","_cell_guid":"8ffe76a4-7e56-481b-8406-22520b573edc","trusted":false},"cell_type":"code","source":"fft_all = []\nnames = []\nfor direct in dirs:\n    waves = [f for f in os.listdir(join(train_audio_path, direct)) if f.endswith('.wav')]\n    for wav in waves:\n        sample_rate, samples = wavfile.read(train_audio_path + direct + '/' + wav)\n        if samples.shape[0] != sample_rate:\n            samples = np.append(samples, np.zeros((sample_rate - samples.shape[0], )))\n        x, val = custom_fft(samples, sample_rate)\n        fft_all.append(val)\n        names.append(direct + '/' + wav)\n\nfft_all = np.array(fft_all)\n\n# Normalization\nfft_all = (fft_all - np.mean(fft_all, axis=0)) / np.std(fft_all, axis=0)\n\n# Dim reduction\npca = PCA(n_components=3)\nfft_all = pca.fit_transform(fft_all)\n\ndef interactive_3d_plot(data, names):\n    scatt = go.Scatter3d(x=data[:, 0], y=data[:, 1], z=data[:, 2], mode='markers', text=names)\n    data = go.Data([scatt])\n    layout = go.Layout(title=\"Anomaly detection\")\n    figure = go.Figure(data=data, layout=layout)\n    py.iplot(figure)\n    \ninteractive_3d_plot(fft_all, names)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f0807b3e37e3ef955abb1995bfa8cccf2d979fd2","_cell_guid":"7c239d25-c761-4c54-87ed-c47c3f641fdd"},"cell_type":"markdown","source":"Notice that there are *yes/e4b02540_nohash_0.wav*, *go/0487ba9b_nohash_0.wav* and more points, that lie far away from the rest. Let's listen to them."},{"metadata":{"collapsed":true,"_uuid":"8f2b077bda2f985ee5dabfa66778d1c3001cb24b","_cell_guid":"70c6eba5-0621-4292-b9bd-b92bc5066ce3","trusted":false},"cell_type":"code","source":"print('Recording go/0487ba9b_nohash_0.wav')\nipd.Audio(join(train_audio_path, 'go/0487ba9b_nohash_0.wav'))","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_uuid":"4116f932a79f3e2e9f745ec4e441c0529ccf66fa","_cell_guid":"1be3b686-d99d-4a1b-aa23-165ad4e5c67e","trusted":false},"cell_type":"code","source":"print('Recording yes/e4b02540_nohash_0.wav')\nipd.Audio(join(train_audio_path, 'yes/e4b02540_nohash_0.wav'))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f1252a107c341e3592a301398be2a56202b9f2a5","_cell_guid":"755472e7-6c19-4c2b-bb8a-d161a501c80a"},"cell_type":"markdown","source":"If you will look for anomalies for individual words, you can find for example this file for *seven*:"},{"metadata":{"collapsed":true,"_uuid":"542a9b0b0c61e2420af77d0f0196178475952681","_cell_guid":"45552683-e626-4143-88c2-b22f9fa00737","trusted":false},"cell_type":"code","source":"print('Recording seven/e4b02540_nohash_0.wav')\nipd.Audio(join(train_audio_path, 'seven/b1114e4f_nohash_0.wav'))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8e532f1e54b04ba5850a163428d1fd94c47bee37","_cell_guid":"67b4915a-7459-4863-a696-ac8220da5a29"},"cell_type":"markdown","source":"That's nothing obviously important. Usually you can find some distortions using this method. Data seems to contain what it should."},{"metadata":{"_uuid":"a50ecff1ab0f1629bbd069c5e325e09035bc2778","_cell_guid":"8e1ba6de-b802-417a-b98d-72d0fed93296"},"cell_type":"markdown","source":"## 3. Where to look for the inspiration\n<a id=\"wheretostart\"></a> \n\nYou can take many different approches for the competition. I can't really advice any of that. I'd like to share my initial thoughts.\n\nThere is a trend in recent years to propose solutions based on neural networks. Usually there are two architectures. My ideas are here.\n\n1. Encoder-decoder: https://arxiv.org/abs/1508.01211\n2. RNNs with CTC loss: https://arxiv.org/abs/1412.5567<br>\nFor me, 1 and 2  are a sensible choice for this competition, especially if you do not have background in SR field. They try to be end-to-end solutions. Speech recognition is a really big topic and it would be hard to get to know important things in short time.\n\n3. Classic speech recognition is described here: http://www.ece.ucsb.edu/Faculty/Rabiner/ece259/Reprints/tutorial%20on%20hmm%20and%20applications.pdf\n\nYou can find *Kaldi* [Tutorial for dummies](http://kaldi-asr.org/doc/kaldi_for_dummies.html), with a problem similar to this competition in some way.\n\n4. Very deep CNN - Don't know if it is used for SR. However, most papers concern Large Vocabulary Continuous Speech Recognition Systems (LVCSR). We got different task here - a very small vocabulary, and recordings with only one word in it, with a (mostly) given length. I suppose such approach can win the competition. \n"},{"metadata":{"_uuid":"f4862ca5ef5f190c0593ff4bd165acc74588c373","_cell_guid":"9192984a-7307-41e5-bb02-9f46f14822c7"},"cell_type":"markdown","source":"---\n\n**If you like my work please upvote.**\n\nLeave a feedback that will let me improve! "}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","file_extension":".py","codemirror_mode":{"version":3,"name":"ipython"},"version":"3.6.3","nbconvert_exporter":"python","mimetype":"text/x-python","pygments_lexer":"ipython3"}},"nbformat":4,"nbformat_minor":1}