{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Wavelets\nEveryone knows the Fourier transform where (stationary) signals are represented as a sum of sinusoids, and the standard transform is only localized in frequency. [Wavelets](https://en.wikipedia.org/wiki/Wavelet) on the other hand are localized in both time and frequency. Wavelets have proven to perform well in many ML applications on non-stationary signals. The two main approaches of using wavelets are:\n* Discrete wavelet transform (DWT): A signal is broken down into subbands, from which features are extracted and used as input to ML models.\n* Continuous wavelet transform (CWT): A scaleogram can be constructed from the signal. The scaleogram can be compared to a spectrogram when using FFT. Scaleogram images can then be used as input to an image classification model.  \n\nIn this notebook we will look into how to create scaleograms from bird song signatures.","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport cv2\nimport matplotlib.image as mpimg\nimport io\nimport librosa\nimport warnings\nwarnings.filterwarnings('ignore') # get rid of librosa warnings","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We could use [PyWavelets](https://pywavelets.readthedocs.io/en/latest/) for the CWT, but here we will use a modified version of [Alexander Neergaard's CWT](https://github.com/aneergaard/CWT). Specifically a parameter is added to allow cutting off the longest wavelets (lowest frequencies).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import math\n\n# Continous Wavelet Transform with Morlet wavelet \n# Original code by Alexander Neergaard, https://github.com/neergaard/CWT\n# \n# Parameters:\n#   data: input data\n#   nv: # of voices (scales) per octave\n#   sr: sampling frequency (Hz)\n#   low_freq: lowest frequency (Hz) of interest (limts longest scale)\ndef cwt2(data, nv=10, sr=1., low_freq=0.):\n    data -= np.mean(data)\n    n_orig = data.size\n    ds = 1 / nv\n    dt = 1 / sr\n\n    # Pad data symmetrically\n    padvalue = n_orig // 2\n    x = np.concatenate((np.flipud(data[0:padvalue]), data, np.flipud(data[-padvalue:])))\n    n = x.size\n\n    # Define scales\n    _, _, wavscales = getDefaultScales(n_orig, ds, sr, low_freq)\n    num_scales = wavscales.size\n\n    # Frequency vector sampling the Fourier transform of the wavelet\n    omega = np.arange(1, math.floor(n / 2) + 1, dtype=np.float64)\n    omega *= (2 * np.pi) / n\n    omega = np.concatenate((np.array([0]), omega, -omega[np.arange(math.floor((n - 1) / 2), 0, -1, dtype=int) - 1]))\n\n    # Compute FFT of the (padded) time series\n    f = np.fft.fft(x)\n\n    # Loop through all the scales and compute wavelet Fourier transform\n    psift, freq = waveft(omega, wavscales)\n\n    # Inverse transform to obtain the wavelet coefficients.\n    cwtcfs = np.fft.ifft(np.kron(np.ones([num_scales, 1]), f) * psift)\n    cfs = cwtcfs[:, padvalue:padvalue + n_orig]\n    freq = freq * sr\n\n    return cfs, freq\n\ndef getDefaultScales(n, ds, sr, low_freq):\n    nv = 1 / ds\n    # Smallest useful scale (default 2 for Morlet)\n    s0 = 2\n\n    # Determine longest useful scale for wavelet\n    max_scale = n // (np.sqrt(2) * s0)\n    if max_scale <= 1:\n        max_scale = n // 2\n    max_scale = np.floor(nv * np.log2(max_scale)) \n    a0 = 2 ** ds\n    scales = s0 * a0 ** np.arange(0, max_scale + 1)\n    \n    # filter out scales below low_freq\n    fourier_factor = 6 / (2 * np.pi)\n    frequencies = sr * fourier_factor / scales\n    frequencies = frequencies[frequencies >= low_freq]\n    scales = scales[0:len(frequencies)]\n\n    return s0, ds, scales\n\ndef waveft(omega, scales):\n    num_freq = omega.size\n    num_scales = scales.size\n    wft = np.zeros([num_scales, num_freq])\n\n    gC = 6\n    mul = 2\n    for jj, scale in enumerate(scales):\n        expnt = -(scale * omega - gC) ** 2 / 2 * (omega > 0)\n        wft[jj, ] = mul * np.exp(expnt) * (omega > 0)\n\n    fourier_factor = gC / (2 * np.pi)\n    frequencies = fourier_factor / scales\n\n    return wft, frequencies","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## FFT vs Spectrogram vs Scaleogram\nFirst we will take a look at the difference between the standard FFT, the spectrogram and the scaleogram.","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import scipy.fftpack\nfrom scipy import signal\nfrom scipy.signal import chirp\n\nDB_RANGE = 100 # dynamic range to show in dB\nSR = 22050\nCMAP = 'magma'\n\ndef show_sigx3(d):\n    fig, axes = plt.subplots(1, 3, figsize=(16,5))\n    # FFT\n    N, T = SR, 1./SR\n    x = np.linspace(0.0, int(N*T), N)\n    yf = scipy.fftpack.fft(d*np.hamming(len(d)))\n    xf = np.linspace(0.0, 1.0/(2.0*T), N//2)\n    axes[0].plot(xf,  20*np.log10(2.0/N *np.abs(yf[:N//2])))\n    axes[0].set_ylim(-80,0)\n    axes[0].set_title('FFT')\n    # Spectrogram\n    f, t, Sxx = signal.spectrogram(d, SR)\n    axes[1].pcolormesh(t, f, 20*np.log10(Sxx), shading='auto', cmap=CMAP, vmax=-60, vmin=-60-DB_RANGE)\n    axes[1].set_title('Spectrogram')\n    # CWT\n    #cs, f = calc_cwt(d)\n    cs, f = cwt2(d, nv=12, sr=SR, low_freq=40)\n    axes[2].imshow(20*np.log10(np.abs(cs)), cmap=CMAP, aspect='auto', norm=None, vmax=0, vmin=-DB_RANGE)\n    axes[2].set_title('Scaleogram')\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Random noise","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"d = np.random.randn(SR)\nshow_sigx3(d)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The scaleogram has high temporal resolution but poor frequency resolution at small scales (high frequencies), and high frequency resolution but poor temporal resolution at large scales (low frequencies).","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Chirp","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"t = np.linspace(0, 1, SR)\nd = chirp(t, f0=50, f1=10e3, t1=1., method='hyperbolic')\nshow_sigx3(d)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Already we can see that scaleograms generally appears smoother than spectrograms, and with less artifacts. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Bird songs\n\nRight, let's look at some bird song recordings. Will resample all recordings to mono, 22050Hz.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from IPython.display import Audio\nfrom sklearn.preprocessing import minmax_scale\n\nTRAIN_DIR = '../input/birdsong-recognition/train_audio'\nSR = 22050 # sample rate in Hz\n\ndef rd_file(fname, offset=0, duration=60):\n    data, _ = librosa.load(fname, sr=SR, mono=True, offset=offset, duration=duration)\n    data = minmax_scale(data-data.mean(), feature_range=(-1,1))\n    return data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# load a file with Black-throated Blue Warbler\nd1 = rd_file(TRAIN_DIR+'/btbwar/XC139608.mp3')\nAudio(d1, rate=SR)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# load another file with Black-throated Green Warbler\nd2 = rd_file(TRAIN_DIR+'/btnwar/XC135117.mp3')\nAudio(d2, rate=SR)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Now, make a wavelet transform of the audio recordings and plot the scaleograms for these two files. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"cs1, f1 = cwt2(d1, nv=12, sr=SR, low_freq=40)\nplt.figure(figsize = (16,4))\nplt.imshow(20*np.log10(np.abs(cs1)), cmap=CMAP, aspect='auto', norm=None, vmax=0, vmin=-30)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cs2, f2 = cwt2(d2, nv=12, sr=SR, low_freq=40)\nplt.figure(figsize = (16,4))\nplt.imshow(20*np.log10(np.abs(cs2)), cmap=CMAP, aspect='auto', norm=None, vmax=0, vmin=-30)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can clearly see the calls. Let's zoom in on a single call from the two species.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_sigx2(d1, d2, name1='data 1', name2='data 2', cwt=True, db_range=30):\n    fig, axes = plt.subplots(1, 2, figsize=(16,4))\n    d = [d1, d2]\n    name = [name1, name2]\n    for i in range(2):\n        if cwt == True:\n            cs, _ = cwt2(d[i], nv=12, sr=SR, low_freq=40)\n            axes[i].imshow(20*np.log10(np.abs(cs)),  cmap=CMAP, aspect='auto', norm=None, vmax=0, vmin=-db_range)\n        else:\n            f, t, Sxx = signal.spectrogram(d[i], SR)\n            axes[i].pcolormesh(t, f, 20*np.log10(Sxx), shading='auto', cmap=CMAP, vmax=-60, vmin=-60-db_range)\n        axes[i].set_title(name[i])\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_sigx2(d1[25000:65000], d2[45000:85000], \n           name1='Black-throated Blue Warbler', name2='Black-throated Green Warbler')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A clear difference! For comparison, let's look at the spectrograms of the same signals.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_sigx2(d1[25000:65000], d2[45000:85000], \n           name1='Black-throated Blue Warbler', name2='Black-throated Green Warbler', cwt=False, db_range=60)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Notice that we are using a threshold function (coefficients below -30dB are set to zero) on the amplitudes of the coefficients here to remove background noise. If we include a larger dynamic range, plots look quite different:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_sigx2(d1[25000:65000], d2[45000:85000], \n           name1='Black-throated Blue Warbler', name2='Black-throated Green Warbler', db_range=120)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_sigx2(d1[25000:65000], d2[45000:85000], \n           name1='Black-throated Blue Warbler', name2='Black-throated Green Warbler', cwt=False, db_range=120)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"One major difference between the Fourier transform and wavelet transform is that using the latter we do not have to worry about using a window function to avoid discontinuities, since wavelets are not continuous functions. We can illustrate this by loading a 10s recording, and then compare the coefficients we get by doing the transform on one 10s segment with 10 times 1s segments concatenated. Let's try a recording of the Northern Mockingbird - a bird that basically never makes the same song twice...","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"d3 = rd_file(TRAIN_DIR+'/normoc/XC54018.mp3', offset=127, duration=10)\nAudio(d3, rate=SR)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cs3a, f3 = cwt2(d3, nv=12, sr=SR, low_freq=40)\ncs3b=np.zeros((97,10*SR), dtype=np.complex)\nfor i in range(10):\n    cs3b[:,i*SR:(i+1)*SR], _ = cwt2(d3[i*SR:(i+1)*SR], nv=12, sr=SR, low_freq=40)\nfig, axes = plt.subplots(3, 1, figsize=(16,12))\naxes[0].imshow(20*np.log10(np.abs(cs3a)), cmap=CMAP, aspect='auto', norm=None, vmax=0, vmin=-40)\naxes[0].set_title('Scaleogram - 10s')\naxes[1].imshow(20*np.log10(np.abs(cs3b)), cmap=CMAP, aspect='auto', norm=None, vmax=0, vmin=-40)\naxes[1].set_title('Scaleogram - 10 x 1s')\ndiff = np.abs(np.abs(cs3a)-np.abs(cs3b))\ndiff[diff == 0.] = 1e-15\naxes[2].imshow(20*np.log10(diff), cmap=CMAP, aspect='auto', norm=None, vmax=0, vmin=-40)\naxes[2].set_title('Scaleogram - difference')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The scaleograms are identical - we can apply the CWT on a signal length optimized for available CPU power and memory, and split and concatenate the coefficients whichever way we want! Now, repeat same exercise with spectrogram.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"f, t, Sxxa = signal.spectrogram(d3, SR)\nSxxb=np.zeros((129,10*98))\nfor i in range(10):\n    _, _, Sxxb[:,i*98:(i+1)*98] = signal.spectrogram(d3[i*SR:(i+1)*SR])\nfig, axes = plt.subplots(3, 1, figsize=(16,12))\naxes[0].pcolormesh(t[0:980], f, 20*np.log10(Sxxa[:,0:980]), shading='auto', cmap=CMAP, vmax=-90, vmin=-130)\naxes[0].set_title('Spectrogram - 10s')\naxes[1].pcolormesh(t[0:980], f, 20*np.log10(Sxxb), shading='auto', cmap=CMAP, vmax=-0, vmin=-40)\naxes[1].set_title('Spectrogram - 10 x 1s')\ndiff = np.abs(np.abs(Sxxa[:,0:980]*3e4)-np.abs(Sxxb))\ndiff[diff == 0.] = 1e-15\naxes[2].pcolormesh(t[0:980], f, 20*np.log10(diff), shading='auto', cmap=CMAP, vmax=-0, vmin=-40)\naxes[2].set_title('Spectrogram - difference')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As, we can see, the spectrogram is a lot more complicated - coefficient scaling is different and we do not even get the same number of coefficients for a 10s as with 10x 1s. A summary so far:  \n* Spectrogram (based on FFT): Easy to understand, hard to use\n* Scaleogram (based on CWT): Hard to understand, easy to use.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"So, what we need now is an automated way of extracting these bird signatures, and dump the signatures to jpeg files as input to an image classification model. There are many ways to locate signatures in recordings, and a basic, rule-based one is used below. Of course, these are easy birds, and more difficult birds (or noisy recordings) present a whole range of challenges...","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Automatic signature extraction\nIf we calculate the variance of each wavelet vector, we might use that as input to the extraction algorithm.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# calculate variance of coefficients\ndef calc_var(cs, thres):\n    c = 20*np.log10(np.abs(cs))\n    c[c < thres] = 0.\n    e = np.var(c, axis=0)\n    return e / max(e)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's see how that works out for the two files.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(2, 1, figsize=(16,4))\nv1 = calc_var(cs1, -30)\naxes[0].plot(v1)\nv2 = calc_var(cs2, -30)\naxes[1].plot(v2)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Scipy has a function find_peaks() that we could use to find high-energy parts of the scaleogram.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.signal import find_peaks\n\ndef mask_sig(n, peaks, sr=22050, dur=0.1):\n    mask = np.zeros(n)\n    subm = int(sr*dur*0.5)\n    if len(peaks > 0):\n        for i in range(len(peaks)):\n            mask[max(peaks[i]-subm, 0): min(peaks[i]+subm, n)] = 1\n    return mask\n\nfig, axes = plt.subplots(2, 1, figsize=(16,4))\n# peak detection + gliding window\npeaks, _ = find_peaks(v1, prominence=0.3)\naxes[0].plot(v1)\naxes[0].plot(peaks, v1[peaks], \"x\")\nm = mask_sig(len(v1), peaks, SR, 0.3)\naxes[1].plot(m)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, we extract the mask segments with other Scipy functions:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_mask(vdata, prom=0.2, dur=0.2, sr=22050):\n    peaks, _ = find_peaks(vdata, prominence=prom)\n    return mask_sig(len(vdata), peaks, sr, dur)\n    \ndef get_regions(mask, sr, species, filename):\n    regions = scipy.ndimage.find_objects(scipy.ndimage.label(mask)[0])\n    regs = []\n    for r in regions:\n        dur = round((r[0].stop-r[0].start)/sr,3)\n        regs.append([r[0].start, r[0].stop, dur,species,filename])\n    return pd.DataFrame(regs, columns=['Start', 'End', 'Duration','Species','File'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mask = get_mask(v1, prom=0.3, dur=0.2, sr=SR)\ndf = get_regions(mask, SR, 'btbwar', 'XC139608.mp3')\ndf.head(6)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We are only interested in segments of a certain length, say 1s or more:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df = df[df.Duration >= 1.0]\ndf = df.reset_index(drop=True)\ndf","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Plot those regions to check that we captured the signatures:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# we plot two seconds\nplot_sigx2(d1[df.Start[0]:df.Start[0]+2*SR], d1[df.Start[1]:df.Start[1]+2*SR], \n           name1='btbwar 1', name2='btbwar 2')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Not bad at all! Let's try the other file:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"mask = get_mask(v2, prom=0.3, dur=0.3, sr=SR)\ndf = get_regions(mask, SR, 'btnwar', 'XC135117.mp3')\ndf = df[df.Duration >= 1.0]\ndf = df.reset_index(drop=True)\ndf","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_sigx2(d2[df.Start[0]:df.Start[0]+2*SR], d2[df.Start[1]:df.Start[1]+2*SR], \n           name1='btnwar 1', name2='btnwar 2')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The second signature is clipped a bit, but that is OK. All five signatures were detected. The final step is to combine everything into a function that will extract signatures from a file and dump them as jpeg images.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"A rule-based signature extraction algorithm is pretty much doomed on complex or noisy signals. A better solution is to train a binary classifier to extract signatures on a reasonable scale (1-5s), and then feed the extracted signatures to a species classifier. Breaking this problem down into two steps simplifies things a lot, and each step can be optimized independently. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# free up memory\n%xdel d1\n%xdel d2\n%xdel cs1\n%xdel cs2","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Create images","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def img_resize(cs, w=512, h=512, log=True, lthres=-30):\n    buf = io.BytesIO()\n    if log == True:\n        plt.imsave(buf, 20*np.log10(np.abs(cs)), cmap=CMAP, format='png', vmax=0, vmin=lthres)\n    else:\n        plt.imsave(buf, np.abs(cs), cmap=CMAP, format='png')\n    buf.seek(0)\n    img_bytes = np.asarray(bytearray(buf.read()), dtype=np.uint8)\n    img = cv2.imdecode(img_bytes, cv2.IMREAD_COLOR)\n    return cv2.resize(img, (w, h), interpolation=cv2.INTER_NEAREST)\n\n# Parameters:\n#   filename: mp3-file\n#   voices: # of scales per octave\n#   sr: sampling frequency (Hz)\n#   low_freq: low freq cutoff (Hz)\n#   thres: scaleogram threshold (dB)\n#   prom: peak detect prominence\n#   peakdur: peak extension (s)\n#   sigthres: smallest signature detection to process (s)\n#   siglen: length of output signature (s)\n#   img_size: output image size\n#   outdir: output directory\ndef scaleo_extract(filename, voices=12, sr=22050, low_freq=40, thres=-30, prom=0.3, \n                   peakdur=0.3, sigthres=1, siglen=2, img_size=512, outdir='.'):\n    d = rd_file(filename)\n    cs, _ = cwt2(d, nv=voices, sr=sr, low_freq=low_freq) # wavelet transform\n    v = calc_var(cs, thres) # coefficient variance\n    peaks, _ = find_peaks(v, prominence=prom) \n    m = mask_sig(len(v), peaks, sr=sr, dur=peakdur) # create signal mask\n    df = get_regions(m, sr, filename.split('/')[-2], filename.split('/')[-1])\n    df = df[df.Duration >= sigthres] # filter out insignificant signatures\n    df = df.reset_index(drop=True)\n    if len(df) > 0:\n        for i in range(len(df)):\n            img = img_resize(cs[:,df.Start[i]:df.Start[i]+siglen*sr], \n                             w=img_size, h=img_size, log=True, lthres=thres)\n            fn = df.Species[i]+'-'+filename.split('/')[-1].split('.')[-2]+\"-{:03d}.jpg\".format(i) \n            cv2.imwrite(outdir+'/'+fn, img)\n    return df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Process a few files and see what we get.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"%mkdir ./tmp","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"flist = ['/btbwar/XC139608.mp3', '/btbwar/XC51863.mp3', '/btbwar/XC134502.mp3', '/btbwar/XC415596.mp3']\nfor i in flist:\n    scaleo_extract(TRAIN_DIR+i, outdir='./tmp')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import glob\n\nimages = glob.glob('./tmp' + \"/btbwar*.jpg\")\nplt.figure(figsize=(20,20))\ncolumns = 4\nfor i, image in enumerate(images):\n    plt.subplot(len(images) / columns + 1, columns, i + 1)\n    plt.imshow(mpimg.imread(image))\n    plt.axis('off')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, another warbler.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"flist = ['/btnwar/XC135117.mp3', '/btnwar/XC135495.mp3', '/btnwar/XC173264.mp3', '/btnwar/XC501261.mp3', '/btnwar/XC486860.mp3']\nfor i in flist:\n    scaleo_extract(TRAIN_DIR+i, outdir='./tmp')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"images = glob.glob('./tmp' + \"/btnwar*.jpg\")\nplt.figure(figsize=(20,25))\ncolumns = 5\nfor i, image in enumerate(images):\n    plt.subplot(len(images) / columns + 1, columns, i + 1)\n    plt.imshow(mpimg.imread(image))\n    plt.axis('off')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We extracted 4 non-signatures from file XC501261 - should not be a problem to do a post check to detect and remove false signatures. Or better: Train an ML model (binary classifier) to detect signatures.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Color or grayscale?\nThe scaleograms represent magnitude of the wavelet coefficients. Pseudo-coloring is used here for human perception. But there is really not more information here than in a grayscale image. Which colormap to use for training? Well, if we think of augmentation with brightness or contrast adjustments, those are only \"correct\" for grayscale. On the other hand, image models might train better with color inputs... So maybe we need to try both - or even an ensemble of the two!","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Summary\nThis notebook shows how to use the contiunous wavelet transform to create scaleograms from non-stationary audio signals, with the purpose of using an image classifier to detect bird species. The outline of a two-stage ML pipeline is starting to emerge:  \n* Step 1: Train a binary classifier to detect and extract bird songs on a reasonable scale (1-5s)\n* Step 2: Feed the bird song segments to a species classifier  \n\nAll using the wavelet transform and scaleograms of course.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"!rm -fr tmp","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}