Cora Johnson-Roberson

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

from skimage import feature, filters, io, morphology
from scipy import ndimage
import librosa
import librosa.display # must be imported separately

from IPython.display import Audio, Image

# Load the Boston skyline image https://commons.wikimedia.org/wiki/File:Boston_Twilight_Panorama_3.jpg
# by Wikimedia Commons user Fcb981 https://commons.wikimedia.org/wiki/User:Fcb981
# CC BY-SA 3.0

skyline_url = "https://upload.wikimedia.org/wikipedia/commons/6/67/Boston_Twilight_Panorama_3.jpg"
skyline = io.imread(skyline_url)
Image(url=skyline_url)

# Cropping and thresholding to capture just the outline.

lx, ly, _ = skyline.shape
skyline_cropped = skyline[: -lx // 4, :, 2]

threshold = filters.threshold_isodata(skyline_cropped)
thresholded = skyline_cropped > threshold

outline = np.argmin(thresholded, axis=0)
outline = skyline_cropped.shape[0] - outline

thresholded = thresholded[::-1, :]

fig, axes = plt.subplots()
axes.plot(outline, color='r')
axes.imshow(thresholded, cmap=plt.cm.gray, origin='lower')

envelope = outline / outline.max()
stretched_envelope = envelope.repeat(10)

# generate white noise in the range [-1.0, 1.0)
waveform = ((np.random.rand(*stretched_envelope.shape) - 0.5) * 2)

# impose envelope on noise
waveform *= stretched_envelope

fig, axes = plt.subplots()
axes.plot(waveform)

Audio(waveform, rate=22050)

fig, axes = plt.subplots()
axes.plot(stretched_envelope)

Audio(stretched_envelope, rate=22050)

skyline_ifft = librosa.istft(thresholded)

D = librosa.stft(skyline_ifft)

fig, axes = plt.subplots()
plt.sca(axes)

librosa.display.specshow(librosa.amplitude_to_db(D, ref=np.max),
                         y_axis='log', x_axis='time')
plt.title('Power spectrogram')
plt.colorbar(format='%+2.0f dB')
plt.tight_layout()

Audio(skyline_ifft, rate=22050)

outline_2d = np.zeros(thresholded.shape)

for i, y in enumerate(outline):
    outline_2d[y, i] = 1.0

outline_ifft = librosa.istft(outline_2d)

D = librosa.stft(outline_ifft)

fig, axes = plt.subplots()
plt.sca(axes)

librosa.display.specshow(librosa.amplitude_to_db(D, ref=np.max),
                         y_axis='log', x_axis='time')
plt.title('Power spectrogram')
plt.colorbar(format='%+2.0f dB')
plt.tight_layout()

Audio(outline_ifft, rate=22050)

import tables

h5file = tables.openFile("msd_summary_file.h5")

h5file

musicbrainz = h5file.root.musicbrainz.songs
analysis = h5file.root.analysis.songs

from collections import defaultdict
by_year = defaultdict(list)
for row1 in musicbrainz.where('year != 0'):
    row2 = analysis[row1.nrow]
    by_year[row1['year']].append(row2['tempo'])
    

print("Year\tAvg. BPM")
for year, values in by_year.iteritems():
    num_songs = len(values)
    print("{}\t{:0.2f} ({} songs)".format(year, float(sum(values))/num_songs, num_songs))

# lowess function from https://gist.github.com/agramfort/850437
import numpy as np
from math import ceil
from scipy import linalg

def lowess(x, y, f=2./3., iter=3):
    """lowess(x, y, f=2./3., iter=3) -> yest

    Lowess smoother: Robust locally weighted regression.
    The lowess function fits a nonparametric regression curve to a scatterplot.
    The arrays x and y contain an equal number of elements; each pair
    (x[i], y[i]) defines a data point in the scatterplot. The function returns
    the estimated (smooth) values of y.

    The smoothing span is given by f. A larger value for f will result in a
    smoother curve. The number of robustifying iterations is given by iter. The
    function will run faster with a smaller number of iterations."""
    n = len(x)
    r = int(ceil(f*n))
    h = [np.sort(np.abs(x - x[i]))[r] for i in range(n)]
    w = np.clip(np.abs((x[:,None] - x[None,:]) / h), 0.0, 1.0)
    w = (1 - w**3)**3
    yest = np.zeros(n)
    delta = np.ones(n)
    for iteration in range(iter):
        for i in range(n):
            weights = delta * w[:,i]
            b = np.array([np.sum(weights*y), np.sum(weights*y*x)])
            A = np.array([[np.sum(weights), np.sum(weights*x)],
                   [np.sum(weights*x), np.sum(weights*x*x)]])
            beta = linalg.solve(A, b)
            yest[i] = beta[0] + beta[1]*x[i]

        residuals = y - yest
        s = np.median(np.abs(residuals))
        delta = np.clip(residuals / (6.0 * s), -1, 1)
        delta = (1 - delta**2)**2

    return yest

%matplotlib notebook
import matplotlib.pyplot as plt

x = np.arange(1950, 2010)
y = np.array([np.mean(by_year[i]) for i in x])

# need to rescale, else we get a singular matrix
scaled_x = np.arange(0,1, 1.0/len(x))
ymax = float(np.max(y))
scaled_y = y / ymax

f = 0.25 # smoothing span; larger == smoother
yest = lowess(scaled_x, scaled_y, f=f, iter=3)

yest *= ymax

plt.xlim(1950, 2010)
plt.scatter(x, y, label='y raw')
plt.plot(x, yest, label='y smoothed')
plt.xlabel('year')
plt.ylabel('Average BPM')
plt.title('BPMs over time')