{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Birds and Frog species identification through deep learning"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nimport warnings\nimport math\nwarnings.filterwarnings('ignore')\nimport os\nimport pandas as pd\n!pip install librosa\nimport librosa\nimport librosa.display\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport IPython.display as ipd\nfrom IPython.core.display import HTML\nfrom IPython.display import Audio\nfrom scipy.io import wavfile\nfrom IPython.display import Audio\nimport csv\nimport pickle\nimport joblib\nfrom sklearn.preprocessing import normalize\nfrom sklearn import preprocessing\nfrom sklearn.ensemble import RandomForestClassifier \nfrom sklearn.metrics import accuracy_score, confusion_matrix, classification_report\nfrom sklearn import preprocessing\nfrom sklearn.model_selection import train_test_split\nimport tensorflow as tf\nimport tensorflow_datasets as tfds","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"os.chdir('/kaggle/input/rfcx-species-audio-detection/')\ndf = pd.read_csv('train_tp.csv') #Load the true positives","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## The max frequency of the signal is 13,687.5 Hz, so use a sample rate that will preserve 14K Hz or 28,000."},{"metadata":{"trusted":true},"cell_type":"code","source":"df.f_max.max()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Look at a sample to see what it sounds like:"},{"metadata":{"trusted":true},"cell_type":"code","source":"os.chdir('/kaggle/input/rfcx-species-audio-detection/train')\nsample_num = 3 #pick a file index\nfilename=df.recording_id[sample_num]+str('.flac') #get the filename\ntstart = df.t_min[sample_num] #define the beginning time of the signal\ntend = df.t_max[sample_num] #define the end time of the signal\ntstart, tend, df.species_id[sample_num],df.songtype_id[sample_num] #print start and end times, the specied_id & songtype_id\ny, sr = librosa.load(filename,sr=28000) #Load the file\ny.shape #print the size of the audio file","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Cut the sample to where the species is recorded:"},{"metadata":{"trusted":true},"cell_type":"code","source":"y_cut=y[int(round(tstart*sr)):int(round(tend*sr))]\nlibrosa.display.waveplot(y_cut,sr=sr, x_axis='time', offset=0.0)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## A sample sounds like this:"},{"metadata":{"trusted":true},"cell_type":"code","source":"ipd.Audio(y_cut, rate=sr) # load a NumPy array\nwavfile.write('/kaggle/working/filename.wav', sr, y_cut)\nAudio('/kaggle/working/filename.wav')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Features for extraction:"},{"metadata":{},"cell_type":"markdown","source":"## 1) Power spectrum"},{"metadata":{"trusted":true},"cell_type":"code","source":"#FFT\nfft = np.fft.fft(y_cut)\nmagnitude=np.abs(fft)\nfrequency = np.linspace(0, sr, len(magnitude))[:int(len(magnitude)/2)]\nmagnitude=magnitude[:int(len(magnitude)/2)]\n# plot spectrum\nplt.figure(figsize=(20,7))\nplt.plot(frequency, magnitude, alpha=0.4)\nplt.xlabel(\"Frequency\",fontsize=15)\nplt.ylabel(\"Magnitude\",fontsize=15)\nplt.title(\"Power spectrum\",fontsize=20)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 2) Spectrogram"},{"metadata":{"trusted":true},"cell_type":"code","source":"hop_length = 512 # in num. of samples\nn_fft = 255 # window in num. of samples\n# perform stft\nstft = librosa.stft(y_cut, n_fft=n_fft, hop_length=hop_length)\n# calculate abs values on complex numbers to get magnitude\nspectrogram = np.abs(stft)\n\n# display spectrogram\nplt.figure(figsize=(20,7))\nlibrosa.display.specshow(librosa.amplitude_to_db(spectrogram),sr=sr, hop_length=hop_length)\nplt.xlabel(\"Time\", fontsize=15)\nplt.xticks()\nplt.ylabel(\"Frequency\", fontsize=15)\nplt.colorbar()\nplt.title(\"Spectrogram\", fontsize=20)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3) MFCC's"},{"metadata":{"trusted":true},"cell_type":"code","source":"MFCCs = librosa.feature.mfcc(y_cut, n_fft=n_fft, hop_length=hop_length,n_mfcc=128)\nfig, ax = plt.subplots(figsize=(20,7))\nlibrosa.display.specshow(MFCCs,sr=sr, hop_length=hop_length)\nax.set_xlabel('Time', fontsize=15)\nax.set_title('MFCC', size=20)\nplt.colorbar()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 4) Spectral centroid"},{"metadata":{"trusted":true},"cell_type":"code","source":"spec_centroid = librosa.feature.spectral_centroid(y=y_cut, sr=sr)\nfig, ax = plt.subplots(figsize=(15,5))\ntime = librosa.times_like(spec_centroid)\nplt.plot(time,spec_centroid.T)\nax.set_xlabel('Time',fontsize=15)\nax.set_title('Spectral Centroid',size=20)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 4) Chroma temperature:"},{"metadata":{"trusted":true},"cell_type":"code","source":"chroma_stft = librosa.feature.chroma_stft(y=y_cut, sr=sr)\nfig, ax = plt.subplots(figsize=(20,7))\nimg = librosa.display.specshow(chroma_stft, y_axis='chroma', x_axis='time', ax=ax)\nax.set_xlabel('Time',fontsize=15)\nax.set_ylabel('Pitch Class', fontsize=15)\nax.set_title('Chromagram', size=20)\nfig.colorbar(img, ax=ax)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 5) Spectral bandwidth"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(20,7))\nspec_bw = librosa.feature.spectral_bandwidth(y=y_cut, sr=sr)\ntimes = librosa.times_like(spec_bw)\nax.semilogy(times, spec_bw[0], label='Spectral bandwidth')\nax.set_ylabel('Hz',fontsize=15)\nax.set_xlabel('Time', fontsize=15)\nax.set_title('Spectral Bandwidth', size=20)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 6) Spectral rolloff"},{"metadata":{"trusted":true},"cell_type":"code","source":"rolloff = librosa.feature.spectral_rolloff(y=y, sr=sr)\ntimes = librosa.times_like(rolloff)\nfig, ax = plt.subplots(figsize=(20,7))\nplt.plot(times,rolloff[0])\nax.set_xlabel('Time', fontsize=15)\nax.set_title('Spectral Rolloff', fontsize=20)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Shape the data into 4D structures suitable for CNN"},{"metadata":{"trusted":true},"cell_type":"code","source":"image=normalize(np.array([spec_bw]).reshape(1,207))\nimage = normalize(np.append(image,spec_centroid.reshape(1,207), axis=0))\nfor i in range(0,9):\n        image = normalize(np.append(image, spec_bw.reshape(1,207), axis=0))\n        image = normalize(np.append(image, spec_centroid.reshape(1,207), axis=0))\n        image = normalize(np.append(image,chroma_stft.reshape(12,207), axis=0))\nimage.shape\nfig, ax = plt.subplots(figsize=(22,10))\nplt.imshow(image)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image=np.dstack((image,spectrogram.reshape(128,207)))\nimage=np.dstack((image,MFCCs.reshape(128,207)))\nimage.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, (ax1,ax2,ax3) = plt.subplots(1,3, figsize=(25,10))\nax1.imshow(image[:,:,0])\nax2.imshow(image[:,:,1])\nax3.imshow(image[:,:,2])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Some species are undersampled with regard to the minority class:"},{"metadata":{"trusted":true},"cell_type":"code","source":"values = df.species_id.value_counts()\nfig, ax = plt.subplots(figsize=(15,5))\nbarlist=plt.bar(values.index[1:],values[1:], label='Minority samples')\nbarlist=plt.bar(values.index[0],values.iloc[0], label='Majority sample')\nplt.xlabel('Species ID', fontsize=17)\nplt.ylabel('Number of audio recordings', fontsize=17)\nplt.title('Majority vs undersampled species', fontsize=20)\nbarlist[0].set_color('orange')\nplt.legend(fontsize=15)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Augment the undersampled species randomly until enough samples have been generated in three ways: \n## <p> 1) adding noise to the signal, <p>2) including more time in the audio sample and <p>3) repeating the audio signal, <p>"},{"metadata":{},"cell_type":"markdown","source":"## Extract background noise from signal prior to augmentation:"},{"metadata":{"trusted":true},"cell_type":"code","source":"y_noise = y[:int(round(tstart*sr))] #noise before species time starts\nipd.Audio(y_noise[20*sr:24*sr], rate=sr) # the section that just has noise in it\nwavfile.write('/kaggle/working/filename-noise.wav', sr, y_noise[20*sr:24*sr]) #this section is devoid of signal\nAudio('/kaggle/working/filename-noise.wav')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Create a noise dictionary to add to the signal during augmentation:"},{"metadata":{"trusted":true},"cell_type":"code","source":"noise1 = np.concatenate((np.array(y_noise[20*sr:24*sr]),np.array(y_noise[20*sr:24*sr]))) #double it to 8 seconds length\n#Augment the noise by adding random noise\nnoise2= noise1 + (0.005 * np.random.randn(len(noise1)))\n#Subtract random noise to get a third type\nnoise3 = noise1 - (0.005 * np.random.randn(len(noise1)))\n#Flip the noises to get three more:\nnoise4 = np.flip(noise1)\nnoise5 = np.flip(noise2)\nnoise6 = np.flip(noise3)\n#place these noise types into a dictionary\nnoise_dict = {1:noise1,2:noise2,3:noise3,4:noise4,5:noise5,6:noise6}\n\nfig, axs = plt.subplots(nrows=3, ncols=2, sharex=True, sharey=True)\n(ax1, ax2), (ax3,ax4), (ax5,ax6) = axs\nimg1=librosa.display.waveplot(noise1,sr=sr, x_axis='time', ax=ax1,offset=0.0)\nimg2=librosa.display.waveplot(noise2,sr=sr, x_axis='time', offset=0.0, ax=ax2)\nimg3=librosa.display.waveplot(noise3,sr=sr, x_axis='time', offset=0.0, ax=ax3)\nimg4=librosa.display.waveplot(noise4,sr=sr, x_axis='time', ax=ax4,offset=0.0)\nimg5=librosa.display.waveplot(noise5,sr=sr, x_axis='time', offset=0.0, ax=ax5)\nimg6=librosa.display.waveplot(noise6,sr=sr, x_axis='time', offset=0.0, ax=ax6)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Create a function to generate features:"},{"metadata":{"trusted":true},"cell_type":"code","source":"#This code was adapted from Nicolas Gervais on https://stackoverflow.com/questions/59241216/padding-numpy-arrays-to-a-specific-size on 1/10/2021\ndef padding(array, xx, yy):\n    \"\"\"\n    :param array: numpy array\n    :param xx: desired height\n    :param yy: desirex width\n    :return: padded array\n    \"\"\"\n\n    h = array.shape[0]\n    w = array.shape[1]\n\n    a = max((xx - h) // 2,0)\n    aa = max(0,xx - a - h)\n\n    b = max(0,(yy - w) // 2)\n    bb = max(yy - b - w,0)\n\n    return np.pad(array, pad_width=((a, aa), (b, bb)), mode='constant')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def generate_features(y_cut):#generate features for machine learning from signal\n    max_size=1000\n    fft = np.fft.fft(y_cut)\n    mag=np.abs(fft)\n    fft=fft[:int(len(mag)/2)]#save half of fft because it's a duplicate of the other half & leave off DC shift\n    stft = padding(np.abs(librosa.stft(y_cut, n_fft=255, hop_length=hop_length)),128,max_size)\n    MFCCs = padding(librosa.feature.mfcc(y_cut, n_fft=n_fft, hop_length=hop_length,n_mfcc=128),128,max_size)\n    spec_centroid = librosa.feature.spectral_centroid(y=y_cut, sr=sr)\n    chroma_stft = librosa.feature.chroma_stft(y=y_cut, sr=sr)\n    spec_bw = librosa.feature.spectral_bandwidth(y=y_cut, sr=sr)\n    rolloff = librosa.feature.spectral_rolloff(y=y_cut, sr=sr)\n    image=np.array([padding(normalize(spec_bw),1,max_size)]).reshape(1,max_size)\n    image = np.append(image,padding(normalize(spec_centroid),1,max_size), axis=0)\n    for i in range(0,9):\n        image = np.append(image, padding(normalize(spec_bw),1,max_size), axis=0)\n        image = np.append(image, padding(normalize(spec_centroid),1,max_size), axis=0)\n        image = np.append(image, padding(normalize(chroma_stft),12,max_size), axis=0)\n    image=np.dstack((image,np.abs(stft)))\n    image=np.dstack((image,MFCCs))\n    return image","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Output of generate_features are 3D padded images of features "},{"metadata":{"trusted":true},"cell_type":"code","source":"image = generate_features(y_cut)\nfig, (ax1,ax2,ax3) = plt.subplots(1,3, figsize=(25,10))\nax1.imshow(image[:,:,0])\nax2.imshow(image[:,:,1])\nax3.imshow(image[:,:,2])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Create a function to generate minority oversampling for the training data"},{"metadata":{"trusted":true},"cell_type":"code","source":"#Minority oversampling algorithm takes in a df of audio with columns=['recording_id','species_id','songtype_id', \n#'t_min', 'f_min', 't_max','f_max'], and an int MM of the number of desired copies and returns df of augmented features\ndef minority_oversample(MM,df_in): \n    #output = np.zeros((2400,128,1000,3)) #make an empty list to store the features\n    features=[]\n    labels = [] #empty array to store labels\n    #For each species, determine how many augmentations are needed\n    df_in=df_in.reset_index()\n    for i in df_in.species_id.unique():\n        print('species_id:',i)\n        filelist = df_in.loc[df_in.species_id == i].index #all the file indices with the same species_id\n        left_overs = MM % len(filelist) \n        for j in range(0,len(filelist)):\n            num_iter = int(math.floor(MM/len(filelist)))-1\n            loop=num_iter\n            if(left_overs>0): loop = num_iter + 1\n            stop = left_overs - 1\n            filename = df_in.iloc[filelist[j]].recording_id +str('.flac') #get the filename\n            tstart = df_in.iloc[filelist[j]].t_min #define the beginning time of the signal\n            tend = df_in.iloc[filelist[j]].t_max #end of signal\n            recording_id = df_in.iloc[filelist[j]].recording_id \n            species_id = i\n            songtype_id = df_in.iloc[filelist[j]].songtype_id\n            y, sr = librosa.load(filename,sr=28000) #Load the file\n            y_cut=y[int(round(tstart*sr)):int(round(tend*sr))] #cut the file to signal start and end\n            data = generate_features(y_cut) #generate features & output numpy array\n            features.append(data[np.newaxis,...])\n            labels.append(species_id)\n            while(loop>0):\n                num=np.random.randint(low=1, high=18) #generate a random number from 1-18\n                if (num>= 1 and num <= 6):\n                    #If 1-6 are generated, add the corresponding noise from the dictionary to the signal.\n                    temp=noise_dict[num][0:len(y_cut)] + y_cut\n                    data=generate_features(temp) #generate features for noise plus signal\n                    features.append(data[np.newaxis,...])\n                    labels.append(species_id)\n                elif (num >6 and num <= 12): \n                    #If 7-12 are generated, increase the time of the audio to include before the tmin start time of the signal.\n                    estart = max(0,tstart - 1) #start the time one second earlier\n                    temp=y[int(round(estart*sr)):int(round(tend*sr))] \n                    data=generate_features(temp) #generate features for longer audio file\n                    features.append(data[np.newaxis,...])\n                    labels.append(species_id)\n                else:\n                    #If 13-18 are generated, double the signal to have it repeated twice.\n                    temp=np.concatenate([y_cut,y_cut],axis=None) #double the audio signal\n                    data=generate_features(temp)\n                    features.append(data[np.newaxis,...])\n                    labels.append(species_id)\n                num_iter = num_iter - 1\n                loop = loop - 1\n                if num_iter < 0:\n                    left_overs = left_overs - 1\n    output=np.concatenate(features,axis=0)\n    return(np.array(output), labels)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Split the data into train and test"},{"metadata":{"trusted":true},"cell_type":"code","source":"X=df.drop('species_id',axis=1)\ny=df.species_id\n#Split into train and test datasets\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=123, stratify=y)\nX_train.shape, X_test.shape, y_train.shape, y_test.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Oversample the minority training true positives and write to file \n## Note that memory constraints during deep learning prevent from upsampling to an ideal max(y_train.value_counts()) and is instead restricted to 50 samples."},{"metadata":{"trusted":true},"cell_type":"code","source":"#Oversample the true positives training dataset, extract features and write to file\ntrain_features, train_labels = minority_oversample(50,pd.concat([X_train,y_train],axis=1)) #sample to 50 species\nprint(train_features.shape,len(train_labels))\njoblib.dump(train_features, '/kaggle/working/train_features.pkl')\njoblib.dump(train_labels, '/kaggle/working/train_labels.pkl')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Extract features for the test dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_features_noOS(df_in): #no oversampling\n    features=[]\n    labels = [] #empty array to store labels\n    #For each species, determine how many augmentations are needed\n    df_in=df_in.reset_index()\n    for i in df_in.species_id.unique():\n        print('species_id:',i)\n        filelist = df_in.loc[df_in.species_id == i].index #all the file indices with the same species_id\n        for j in range(0,len(filelist)):\n            filename = df_in.iloc[filelist[j]].recording_id +str('.flac') #get the filename\n            tstart = df_in.iloc[filelist[j]].t_min #define the beginning time of the signal\n            tend = df_in.iloc[filelist[j]].t_max #end of signal\n            recording_id = df_in.iloc[filelist[j]].recording_id \n            species_id = i\n            songtype_id = df_in.iloc[filelist[j]].songtype_id\n            y, sr = librosa.load(filename,sr=28000) #Load the file\n            y_cut=y[int(round(tstart*sr)):int(round(tend*sr))] #cut the file to signal start and end\n            data = generate_features(y_cut) #generate features & output numpy array\n            features.append(data[np.newaxis,...])\n            labels.append(species_id)\n    output=np.concatenate(features,axis=0)\n    return(np.array(output), labels)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Extract features for the test set and write to file"},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"test_features, test_labels = get_features_noOS(pd.concat([X_test,y_test],axis=1))\nprint(test_features.shape,len(test_labels))\njoblib.dump(test_features, '/kaggle/working/test_features.pkl')\njoblib.dump(test_labels, '/kaggle/working/test_labels.pkl')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Compare against modeling without augmenting training data"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_noOS_features, train_noOS_labels = get_features_noOS(pd.concat([X_train,y_train],axis=1))\nprint(train_noOS_features.shape,len(train_noOS_labels))\njoblib.dump(train_noOS_features, '/kaggle/working/train_noOS_features.pkl')\njoblib.dump(train_noOS_labels, '/kaggle/working/train_noOS_labels.pkl')","execution_count":null,"outputs":[]}],"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}