{"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":"# Audio -> Spectrogram\n\nLoading a 5-sec clip and converting it to a spectrogram.","metadata":{}},{"cell_type":"code","source":"import librosa as librosa\nimport numpy as np\nimport pandas as pd\nimport glob\nimport matplotlib.pyplot as plt\nimport scipy\n%matplotlib inline\nimport librosa.display\nfrom tqdm.notebook import tqdm\nimport warnings  \nwarnings.filterwarnings('ignore') # Librosa's specshow gives warnings about a matplotlib deprecation thing","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get a list of audio recordings\ntrain_files = glob.glob('/kaggle/input/birdclef-2021/train_short_audio/*/*.ogg')\ntrain_files[0]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We want to go from a 5-second clip to a spectrogram, ideally one that captures as much info as possible. Enter librosa and the mel spectrogram","metadata":{}},{"cell_type":"code","source":"# Load 5 seconds of audio\nchunk, rate = librosa.load(train_files[0], sr=32000, offset=3, duration=5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def chunk_to_spec(chunk, SPEC_HEIGHT=64,SPEC_WIDTH=256, rate=32000, FMIN=500, FMAX=12500):\n    mel_spec = librosa.feature.melspectrogram(y=chunk, \n                                              sr=32000, \n                                              n_fft=1024, \n                                              hop_length=int(32000 * 5 / (SPEC_WIDTH - 1)), \n                                              n_mels=SPEC_HEIGHT, \n                                              fmin=FMIN, \n                                              fmax=FMAX)\n    mel_spec_db = librosa.power_to_db(mel_spec, ref=np.max)\n    return mel_spec_db\n\ndef display_spec(spec, SPEC_HEIGHT=64,SPEC_WIDTH=256, rate=32000, FMIN=500, FMAX=12500):\n    librosa.display.specshow(spec, \n                             sr=32000, \n                             hop_length=int(32000 * 5 / (SPEC_WIDTH - 1)), \n                             x_axis='time', \n                             y_axis='mel',\n                             fmin=FMIN, \n                             fmax=FMAX, \n                             cmap=plt.get_cmap('viridis'))\n    \nspec = chunk_to_spec(chunk)\n# display_spec(spec) # Will give axis labels in Hs and time\nplt.imshow(spec, cmap='inferno') # Simple and compact","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Picking the right chunks for training\n\nWhere in a 2-minute recording are the actual calls? I borrowed a trick from https://github.com/BirdVox/PCEN-SNR, computing a sort of signal-to-noise value over time. We can use this to pick peaks where there is most likely to be a call. You can skip the peak_plot code (it is for visualization) - they key function here is get_peaks","metadata":{}},{"cell_type":"code","source":"def peak_plot(y, sr, FMIN=500, FMAX=12500):\n    # PCEN spec\n    plt.figure(figsize=(10, 6))\n    plt.subplot(3, 1, 1)\n    melspec = librosa.feature.melspectrogram(y, sr=sr,\n        fmin=FMIN, fmax=FMAX, n_mels=64)\n    pcen = librosa.core.pcen(melspec, sr=sr,\n        gain=0.8, bias=10, power=0.25, time_constant=0.06, eps=1e-06)\n    librosa.display.specshow(pcen, sr=sr,\n        fmin=FMIN, fmax=FMAX,\n        x_axis='time', y_axis='mel', cmap='magma_r')\n    # plt.title('PCEN-based SNR')\n    plt.tight_layout()\n\n    # SNR and a smoothed SNR with kernel 15\n    plt.subplot(3, 1, 2)\n    pcen_snr = np.max(pcen,axis=0) - np.min(pcen,axis=0)\n    pcen_snr = librosa.power_to_db(pcen_snr / np.median(pcen_snr))\n    median_pcen_snr = scipy.signal.medfilt(pcen_snr, kernel_size=15)\n    times = np.linspace(0, len(y)/sr, num=melspec.shape[1])\n    plt.plot(times, pcen_snr, color=\"orange\")\n    plt.plot(times, median_pcen_snr, color=\"blue\")\n    plt.xlim(times[0], times[-1])\n    plt.ylim(0, 10)\n    # And go through, picking some peaks\n    for i in range(12):\n        t_peak = np.argmax(median_pcen_snr)\n        plt.scatter(times[t_peak], median_pcen_snr[t_peak], c='red', zorder=100)\n        median_pcen_snr[t_peak-50:t_peak+50] = 0 # zero out around the peak to find the next one\n\n    # Kernel 55 for even smoother (bad idea?)\n    plt.subplot(3, 1, 3)\n    median_pcen_snr = scipy.signal.medfilt(pcen_snr, kernel_size=55)\n    times = np.linspace(0, len(y)/sr, num=melspec.shape[1])\n    plt.plot(times, median_pcen_snr, color=\"blue\")\n    plt.xlim(times[0], times[-1])\n    plt.ylim(0, 10)\n    # And go through, picking some peaks\n    for i in range(12):\n        t_peak = np.argmax(median_pcen_snr)\n        plt.scatter(times[t_peak], median_pcen_snr[t_peak], c='red', zorder=100)\n        median_pcen_snr[t_peak-50:t_peak+50] = 0 # zero out around the peak to find the next one","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y, sr = librosa.load(train_files[0], sr=32000, duration=30) # 30 seconds to play with\npeak_plot(y, sr)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y, sr = librosa.load(train_files[100], sr=32000) # A whole recording\npeak_plot(y, sr)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_peaks(y, sr, n_peaks=12, kernel_size=15, zero_dist=50, FMIN=500, FMAX=12500):\n    # Spec\n    melspec = librosa.feature.melspectrogram(y, sr=sr,\n        fmin=FMIN, fmax=FMAX, n_mels=64)\n    pcen = librosa.core.pcen(melspec, sr=sr,\n        gain=0.8, bias=10, power=0.25, time_constant=0.06, eps=1e-06)\n    # SNR\n    pcen_snr = np.max(pcen,axis=0) - np.min(pcen,axis=0)\n    pcen_snr = librosa.power_to_db(pcen_snr / np.median(pcen_snr))\n    # SMoothed SNR\n    median_pcen_snr = scipy.signal.medfilt(pcen_snr, kernel_size=kernel_size)\n    # And go through, picking some peaks\n    times = np.linspace(0, len(y)/sr, num=melspec.shape[1])\n    peak_locs = []\n    for i in range(n_peaks):\n        t_peak = np.argmax(median_pcen_snr)\n        peak_locs.append(times[t_peak])\n        median_pcen_snr[t_peak-50:t_peak+50] = 0 # zero out around the peak to find the next one\n\n    return peak_locs\n\nprint(get_peaks(y, sr, n_peaks=5))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Storing Peak Locations and Labels for Later\n\nThis takes a while to run, so we save the locations of 20 peaks per recording which we can then use in other notebooks down the road.","metadata":{}},{"cell_type":"code","source":"info = []\nfor f in tqdm(train_files):\n    y, sr = librosa.load(f, sr=32000)\n    peaks =  get_peaks(y, sr, n_peaks=20)\n    info.append({\n        'fn':f,\n        'len':len(y), \n        'label':f.split('train_short_audio/')[1].split('/')[0],\n        'peaks':'#'.join([str(p)[:5] for p in peaks]),  \n    })\n\n# Save to file        \ninfo_df = pd.DataFrame(info)\ninfo_df.to_csv('/kaggle/working/info_df.csv', index=False)\ninfo_df.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This file gets saved in /kaggle/working and can be imported into future notebooks. Nifty :)","metadata":{}},{"cell_type":"code","source":"# Old code: saving a spectrogram for the first peak in each file\n# def getx(item, peakn=0):\n#     f, l, peaks = item['fn'],item['len'], item['peaks']\n#     start_time = min(l/32000 - 5, max(float(peaks.split('#')[peakn])-2.5, 2.5))\n#     y, sr = librosa.load(f, sr=32000, offset=start_time, duration=5)\n#     spec = chunk_to_spec(y, SPEC_HEIGHT=128)\n#     return spec#.astype(np.uint8)\n\n# spec = getx(info_df.sample().to_dict(orient='records')[0])\n# plt.imshow(spec, cmap='inferno')\n\n# !mkdir -p specs\n\n# recs = list(info_df.to_dict(orient='records'))\n# for i in tqdm(range(len(recs))):\n#     row = recs[i]\n#     spec = getx(row, peakn=0)\n#     np.save( f'specs/{i}_0.npy', spec)\n\n# !tar -zcvf specs.tar.gz ./specs > /dev/null\n# !rm -rf specs","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}