{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Birdsong: Simple EDA \n\n\n\n---","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"+ 1. Setup
\n+ 2. Basic EDA
\n+ 3. Basic preprocessing
","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"This competition is used to determine birdsong's bird of origin - it has far reaching implications in studies pertaining to wildlife and conservation. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Setup
","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import IPython.display as ipd # To play sound in the notebook\nfname = '../input/birdsong-recognition/train_audio/aldfly/' + 'XC134874.mp3' # Hi-hat\nipd.Audio(fname)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"!pip install pydub\nfrom pydub import AudioSegment\nimport wave\nfrom scipy.io import wavfile\nimport pandas as pd\nimport re\nimport torch\nimport torchaudio\nfrom torchaudio import transforms\nimport altair as alt\nimport plotly.express as px, plotly.graph_objects as go\nfrom plotly.offline import init_notebook_mode, iplot\ninit_notebook_mode(connected=True)\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\ndef _generate_bar_plot_hor(df, col, title, color, w=None, h=None, lm=0, limit=100):\n cnt_srs = df[col].value_counts()[:limit]\n trace = go.Bar(y=cnt_srs.index[::-1], x=cnt_srs.values[::-1], orientation = 'h',\n marker=dict(color=color))\n\n layout = dict(title=title, margin=dict(l=lm), width=w, height=h)\n data = [trace]\n annotations = []\n annotations += [go.layout.Annotation(x=673, y=100, xref=\"x\", yref=\"y\", text=\"(Most Popular)\", showarrow=False, arrowhead=7, ax=0, ay=-40)]\n fig = go.Figure(data=data, layout=layout)\n iplot(fig)\nplt.style.use('fivethirtyeight')\nsound = AudioSegment.from_mp3(fname)\nsound.export('fex1.wav', format=\"wav\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now that we have our example file in .wav format, we can proceed with our analysis:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"wav = wave.open('fex1.wav')\nrate, data = wavfile.read('fex1.wav')\nprint(\"Sampling (frame) rate = \", rate)\nprint(\"Total samples (frames) = \", data.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can take a look at this with, well, good ol' Matplotlib.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Basic EDA
","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(data, '-', );","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have a few jolts at the end (could be the call finally becoming louder, approacing a *crescendo*).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 4))\nplt.plot(data[:500], '.'); plt.plot(data[:500], '-');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We see that it is fairly inconsistent in first 500 frames, with a few jolts clustered towards the end. This could potientally be a precedent for the true call.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/birdsong-recognition/train.csv')\ntest = pd.read_csv('../input/birdsong-recognition/test.csv')\n\nfig = plt.figure(figsize=(16,5))\ntrain.duration.hist(bins=100);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Majority of files are relatively short it seems.\n\nWe can play some sample audio with IPython's capabilities to display audio in a Jupyter Notebook (very helpful). We will look at the first ten audio files:","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"import IPython\nfor i in train.filename.head(10).values:\n fname = '../input/birdsong-recognition/train_audio/aldfly/' + i\n IPython.display.display(ipd.Audio(fname))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Hear anything? These are the various calls of the Alder Flycatcher - a bird which we are supposed to identify from its call. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"We also have a bunch of other files to examine:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from os import *\ntrain_dir = '../input/birdsong-recognition/train_audio/aldfly/'\ntrain_ = [f for f in listdir(train_dir) if path.isfile(path.join(train_dir, f))]\nfig = plt.figure(figsize=(15,10))\nplt.suptitle('Audio examination among files')\ntrain2_ = []\nfor i in range(1,11):\n fig.add_subplot(5, 2, i)\n fname = train_dir + train_[i]\n sound = AudioSegment.from_mp3(fname)\n train_[i] = re.sub('.mp3', '', train_[i])\n sound.export(f'{i}.wav', format=\"wav\")\n rate, data = wavfile.read(f'{i}.wav')\n plt.plot(data, '-', );\n fig.add_subplot","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The files have incredibly high variance among themselves - the sound is just all over the place. This is very inconsistent and it will prove a challenge for us. We will now move to `duration`:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig,(ax1,ax2)=plt.subplots(1,2,figsize=(15,10))\nax1.hist(train.duration,color='red')\nax1.set_title('duration in train')\nax2.hist(test.seconds,color='blue')\nax2.set_title('duration in test')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Even the durations are inconsistent - this will be quite a task for us. We have lengthy amounts of information in one place but virtually nothing on the other hand. We will try to strip train files of useless information. But first, since we have the goegraphical data, we will plot it on a map.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import folium, numpy as np\nfrom folium.plugins import HeatMap,MarkerCluster\nmy_map_1 = folium.Map(location=[0,0], zoom_start=2, tiles='Stamen Toner')\ntrain2 = train.head(200)\n\nnew_df = pd.DataFrame({\n 'latitude': train2.dropna(axis=0, subset=['latitude','longitude'])['latitude'],\n 'longitude': train2.dropna(axis=0, subset=['latitude','longitude'])['longitude']\n})\nhm = HeatMap(new_df,auto_play=True,max_opacity=0.8)\nhm.add_to(my_map_1)\nmy_map_1 # display","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that the study covers part of South America and most of the USA and Canada. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Basic preprocessing
","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"We will need to trim all the files to a bare minimum length of what is necessary, first we boil it down to 500 frames:","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"from os import *\ntrain_dir = '../input/birdsong-recognition/train_audio/aldfly/'\ntrain_ = [f for f in listdir(train_dir) if path.isfile(path.join(train_dir, f))]\nfig = plt.figure(figsize=(15,10))\nplt.suptitle('Audio examination among files')\ntrain2_ = []\nfor i in range(1,11):\n fig.add_subplot(5, 2, i)\n fname = train_dir + train_[i]\n sound = AudioSegment.from_mp3(fname)\n train_[i] = re.sub('.mp3', '', train_[i])\n sound.export(f'{i}.wav', format=\"wav\")\n rate, data = wavfile.read(f'{i}.wav')\n plt.plot(data[:500], '-', );\n fig.add_subplot","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"No, no, not a good idea. What about the middle 500 frames?","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(15,10))\nplt.suptitle('Audio examination among files')\ntrain2_ = []\nfor i in range(1,11):\n fig.add_subplot(5, 2, i)\n fname = train_dir + train_[i] + '.mp3'\n \n sound = AudioSegment.from_mp3(fname)\n sound.export(f'{i}.wav', format=\"wav\")\n rate, data = wavfile.read(f'{i}.wav')\n plt.plot(data[1000:1500], '-', );\n fig.add_subplot","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"OK, this is better than last 500, but it is not perfect yet. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(15,10))\nplt.suptitle('Audio examination among files')\ntrain2_ = []\nfor i in range(1,11):\n fig.add_subplot(5, 2, i)\n fname = train_dir + train_[i] + '.mp3'\n sound = AudioSegment.from_mp3(fname)\n sound.export(f'{i}.wav', format=\"wav\")\n rate, data = wavfile.read(f'{i}.wav')\n plt.plot(data[1000:2000], '-', );\n fig.add_subplot","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"yes, this is even better, but what about the most inconsistent one? We can use fourier transforms on the audio data:","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"import warnings;warnings.filterwarnings('ignore')\nfrom scipy.fft import fft, ifft\nfig = plt.figure(figsize=(15,10))\nplt.suptitle('Fourier transform among files')\ntrain2_ = []\nfor i in range(1,11):\n fig.add_subplot(5, 2, i)\n fname = train_dir + train_[i] + '.mp3'\n sound = AudioSegment.from_mp3(fname)\n sound.export(f'{i}.wav', format=\"wav\")\n rate, data = wavfile.read(f'{i}.wav')\n plt.plot(fft(data[500:1500]), '-', );\n fig.add_subplot","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In case anyone wants, here's simple yet brilliant Fourier Transform video: https://www.khanacademy.org/science/electrical-engineering/ee-signals/ee-fourier-series/v/ee-fourier-series-intro","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"That is not all! We can use inverse transform as well:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import warnings;warnings.filterwarnings('ignore')\nfrom scipy.fft import fft, ifft\nfig = plt.figure(figsize=(15,10))\nplt.suptitle('Fourier transform among files')\ntrain2_ = []\nfor i in range(1,11):\n fig.add_subplot(5, 2, i)\n fname = train_dir + train_[i] + '.mp3'\n sound = AudioSegment.from_mp3(fname)\n sound.export(f'{i}.wav', format=\"wav\")\n rate, data = wavfile.read(f'{i}.wav')\n plt.plot(ifft(data[500:1500]), '-', );\n fig.add_subplot","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we can \"hose down\" our data with simple signal denoising.","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"import pywt, numpy as np ### FROM https://www.kaggle.com/jackvial/dwt-signal-denoising\ndef maddest(d, axis=None):\n \"\"\"\n Mean Absolute Deviation\n \"\"\"\n return np.mean(np.absolute(d - np.mean(d, axis)), axis)\n\ndef high_pass_filter(x, low_cutoff=1000, sample_rate=0.02):\n \"\"\"\n From @randxie https://github.com/randxie/Kaggle-VSB-Baseline/blob/master/src/utils/util_signal.py\n Modified to work with scipy version 1.1.0 which does not have the fs parameter\n \"\"\"\n \n # nyquist frequency is half the sample rate https://en.wikipedia.org/wiki/Nyquist_frequency\n nyquist = 0.5 * sample_rate\n norm_low_cutoff = low_cutoff / nyquist\n \n # Fault pattern usually exists in high frequency band. According to literature, the pattern is visible above 10^4 Hz.\n # scipy version 1.2.0\n #sos = butter(10, low_freq, btype='hp', fs=sample_fs, output='sos')\n \n # scipy version 1.1.0\n sos = butter(10, Wn=[norm_low_cutoff], btype='highpass', output='sos')\n filtered_sig = signal.sosfilt(sos, x)\n\n return filtered_sig\ndef denoise_signal( x, wavelet='db4', level=1):\n \"\"\"\n 1. Adapted from waveletSmooth function found here:\n http://connor-johnson.com/2016/01/24/using-pywavelets-to-remove-high-frequency-noise/\n 2. Threshold equation and using hard mode in threshold as mentioned\n in section '3.2 denoising based on optimized singular values' from paper by Tomas Vantuch:\n http://dspace.vsb.cz/bitstream/handle/10084/133114/VAN431_FEI_P1807_1801V001_2018.pdf\n \"\"\"\n \n # Decompose to get the wavelet coefficients\n coeff = pywt.wavedec( x, wavelet, mode=\"per\" )\n \n # Calculate sigma for threshold as defined in http://dspace.vsb.cz/bitstream/handle/10084/133114/VAN431_FEI_P1807_1801V001_2018.pdf\n # As noted by @harshit92 MAD referred to in the paper is Mean Absolute Deviation not Median Absolute Deviation\n sigma = (1/0.6745) * maddest( coeff[-level] )\n\n # Calculte the univeral threshold\n uthresh = sigma * np.sqrt( 2*np.log( len( x ) ) )\n coeff[1:] = ( pywt.threshold( i, value=uthresh, mode='hard' ) for i in coeff[1:] )\n \n # Reconstruct the signal using the thresholded coefficients\n return pywt.waverec( coeff, wavelet, mode='per' )\nfig = plt.figure(figsize=(15,10))\nplt.suptitle('Denoise')\ntrain2_ = []\nfor i in range(1,11):\n fig.add_subplot(5, 2, i)\n fname = train_dir + train_[i] + '.mp3'\n sound = AudioSegment.from_mp3(fname)\n sound.export(f'{i}.wav', format=\"wav\")\n rate, data = wavfile.read(f'{i}.wav')\n plt.plot(denoise_signal(data[500:1500]), '-', );\n wavfile.write(f'{i}_lin.wav', 44100,denoise_signal(data[500:1500]))\n fig.add_subplot","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"One with a keen ear can notice the keen differences in the sound which our model should identify. It's your choice: if you want to use it as preprocessing or if you want to use it as data augmentation. Random noise also is possible:","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(15,10))\nplt.suptitle('Denoise')\ntrain2_ = []\nfor i in range(1,11):\n fig.add_subplot(5, 2, i)\n fname = train_dir + train_[i] + '.mp3'\n sound = AudioSegment.from_mp3(fname)\n sound.export(f'{i}.wav', format=\"wav\")\n rate, data = wavfile.read(f'{i}.wav')\n noise = np.random.normal(0, .1, data.shape)\n new = data + noise\n plt.plot(new[500:1500], '-', );\n wavfile.write(f'{i}_lin.wav', 44100,new[500:1500])\n fig.add_subplot","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.DataFrame(train.ebird_code.value_counts())\ndf['name'] = df.index\nimport altair_render_script\nimport altair as alt\nalt.Chart(df).mark_bar().encode(\n x='name',\n y='ebird_code',\n tooltip=[\"name\",\"ebird_code\"]\n).interactive()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A lot of birds feature 100 times in the dataset.","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"df = pd.DataFrame(train.author.value_counts())\ndf['name'] = df.index\nimport altair_render_script\nimport altair as alt\na = alt.Chart(df).mark_bar().encode(\n x='name',\n y='author',\n tooltip=[\"name\",\"author\"]\n).interactive()\ndf = pd.DataFrame(train.ebird_code.value_counts())\ndf['name'] = df.index\nimport altair_render_script\nimport altair as alt\nb = alt.Chart(df).mark_bar().encode(\n x='name',\n y='ebird_code',\n tooltip=[\"name\",\"ebird_code\"]\n).interactive()\na & b","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have a huge amount of authors but one man is responsible for the majority - Paul Marvin.","execution_count":null}],"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":4,"nbformat_minor":4}