rcParams['figure.figsize'] = (16, 4) #wide graphs by default

from scipy.io import wavfile
from scipy.interpolate import interp1d
from scipy.stats import moment
sources = []
sources.append((wavfile.read("media/01 Everlasting Light.wav")))
sources.append((wavfile.read("media/01 P-Funk (Wants to Get Funked Up).wav")))
sources.append((wavfile.read("media/02 fallingElevatorBlues.wav")))
sources.append((wavfile.read("media/02 The Wilhelm Scream.wav")))
sources.append((wavfile.read("media/1-01 Central Texas.wav")))

def windowed_rms(input_sig, win_sizes = [512, 1024, 2048, 4096], hop=None):
    rms_windows = []
    for win_size in win_sizes:
        if not hop:
            hop = win_size/2 
        window_start = arange(0, len(input_sig) - win_size, hop)
        rms = []
        for start in window_start:
            w = input_sig[start: start+win_size].astype(float)
            rms_inst = sqrt(mean(w**2))
            rms.append(rms_inst)
        rms_windows.append(rms)
    return rms_windows, win_sizes

src = sources[0][1]
sr - sources[0][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
in_sig = src_mono
zc = argwhere((in_sig[:-1]*in_sig[1:]) < 0)

zc = r_[zc, argwhere(in_sig == 0)]
plot(in_sig);
plot(zc, zeros_like(zc), 'o');
duration = len(src_mono) / sr
xps = len(zc) / duration
print 'Black Keys - Everlasting Light'
print 'Zero crossings per second: {}'.format(xps)

src = sources[1][1]
sr - sources[1][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
in_sig = src_mono
zc = argwhere((in_sig[:-1]*in_sig[1:]) < 0)

zc = r_[zc, argwhere(in_sig == 0)]
plot(in_sig);
plot(zc, zeros_like(zc), 'o');
duration = len(src_mono) / sr
xps = len(zc) / duration
print 'Parliament - P-Funk (Wants to Get Funked Up)'
print 'Zero crossings per second: {}'.format(xps)

src = sources[2][1]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
#abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
#src_mono_norm = src_mono.astype(float) / abs_max

in_sig = src_mono
zc = argwhere((in_sig[:-1]*in_sig[1:]) < 0)

zc = r_[zc, argwhere(in_sig == 0)]
plot(in_sig);
plot(zc, zeros_like(zc), 'o');
duration = len(src_mono) / sources[0][0]
xps = len(zc) / duration
print 'Me - Falling Elevator Blues'
print 'Zero crossings per second: {}'.format(xps)

src = sources[3][1]
sr - sources[3][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
in_sig = src_mono
zc = argwhere((in_sig[:-1]*in_sig[1:]) < 0)

zc = r_[zc, argwhere(in_sig == 0)]
plot(in_sig);
plot(zc, zeros_like(zc), 'o');
duration = len(src_mono) / sr
xps = len(zc) / duration
print 'James Blake - The Wilhelm Scream'
print 'Zero crossings per second: {}'.format(xps)

src = sources[4][1]
sr - sources[4][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
in_sig = src_mono
zc = argwhere((in_sig[:-1]*in_sig[1:]) < 0)

zc = r_[zc, argwhere(in_sig == 0)]
plot(in_sig);
plot(zc, zeros_like(zc), 'o');
duration = len(src_mono) / sr
xps = len(zc) / duration
print 'Stars of the Lid - Central Texas'
print 'Zero crossings per second: {}'.format(xps)

src = sources[0][1]
sr = sources[0][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max


L = 0.1
win_len = sr * L

w_rms, win_lens = windowed_rms(src_mono_norm, win_sizes=[win_len])

sample_dur_secs = len(src_mono_norm)/float(sr)
time = linspace(0, sample_dur_secs, len(w_rms[0]))

tc = sum(w_rms[0]*time)/sum(w_rms[0])
print tc, "Temporal Centroid"

plot(linspace(0, sample_dur_secs, len(src_mono)), src_mono, alpha=0.5);
vlines(tc, -10000, 10000, lw=4);
grid();

print 'Black Keys - Everlasting Light'
title('Temporal Centroid: {:.2f}% of song'.format(100*(tc/sample_dur_secs)))

src = sources[1][1]
sr = sources[1][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max


L = 0.1
win_len = sr * L

w_rms, win_lens = windowed_rms(src_mono_norm, win_sizes=[win_len])

sample_dur_secs = len(src_mono_norm)/float(sr)
time = linspace(0, sample_dur_secs, len(w_rms[0]))

tc = sum(w_rms[0]*time)/sum(w_rms[0])
print tc, "Temporal Centroid"

plot(linspace(0, sample_dur_secs, len(src_mono)), src_mono, alpha=0.5);
vlines(tc, -10000, 10000, lw=4);
grid();

print 'Parliament - P-Funk (Wants to Get Funked Up)'
title('Temporal Centroid: {:.2f}% of song'.format(100*(tc/sample_dur_secs)))

src = sources[2][1]
sr = sources[2][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max


L = 0.1
win_len = sr * L

w_rms, win_lens = windowed_rms(src_mono_norm, win_sizes=[win_len])

sample_dur_secs = len(src_mono_norm)/float(sr)
time = linspace(0, sample_dur_secs, len(w_rms[0]))

tc = sum(w_rms[0]*time)/sum(w_rms[0])
print tc, "Temporal Centroid"

plot(linspace(0, sample_dur_secs, len(src_mono)), src_mono, alpha=0.5);
vlines(tc, -10000, 10000, lw=4);
grid();

print 'Me - Falling Elevator Blues'
title('Temporal Centroid: {:.2f}% of song'.format(100*(tc/sample_dur_secs)))

src = sources[3][1]
sr = sources[3][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max


L = 0.1
win_len = sr * L

w_rms, win_lens = windowed_rms(src_mono_norm, win_sizes=[win_len])

sample_dur_secs = len(src_mono_norm)/float(sr)
time = linspace(0, sample_dur_secs, len(w_rms[0]))

tc = sum(w_rms[0]*time)/sum(w_rms[0])
print tc, "Temporal Centroid"

plot(linspace(0, sample_dur_secs, len(src_mono)), src_mono, alpha=0.5);
vlines(tc, -10000, 10000, lw=4);
grid();

print 'James Blake - The Wilhelm Scream'
title('Temporal Centroid: {:.2f}% of song'.format(100*(tc/sample_dur_secs)))

src = sources[4][1]
sr = sources[4][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max


L = 0.1
win_len = sr * L

w_rms, win_lens = windowed_rms(src_mono_norm, win_sizes=[win_len])

sample_dur_secs = len(src_mono_norm)/float(sr)
time = linspace(0, sample_dur_secs, len(w_rms[0]))

tc = sum(w_rms[0]*time)/sum(w_rms[0])
print tc, "Temporal Centroid"

plot(linspace(0, sample_dur_secs, len(src_mono)), src_mono, alpha=0.5);
vlines(tc, -10000, 10000, lw=4);
grid();

print 'Stars of the Lid - Central Texas'
title('Temporal Centroid: {:.2f}% of song'.format(100*(tc/sample_dur_secs)))

src = sources[0][1]
sr = sources[0][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=2048, Fs=sr, window=window_hanning, noverlap=1024)
X = sqrt(Pxx)
centroid = []

for spec in X.T:
    sc = sum(spec*freqs)/sum(spec)
    centroid.append(sc)

win_len = 0.1 * sr
w_rms, win_lens = windowed_rms(src_mono_norm, win_sizes=[win_len])

ifunc = interp1d(linspace(0, len(src_mono_norm)/float(sr), len(w_rms[0])), w_rms[0])
interp_rms = ifunc(times)
centroid_rms = array(centroid) * interp_rms
plot(times, centroid_rms, alpha=0.5)
title('Spectral Centroid Scaled by RMS')
mean_sc = sum(centroid) / len(centroid)

print 'Black Keys - Everlasting Light'
print 'Mean Spectral Centroid: {:.2f} Hz'.format(mean_sc)

src = sources[1][1]
sr = sources[1][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=2048, Fs=sr, window=window_hanning, noverlap=1024)
X = sqrt(Pxx)
centroid = []

for spec in X.T:
    sc = sum(spec*freqs)/sum(spec)
    centroid.append(sc)

win_len = 0.1 * sr
w_rms, win_lens = windowed_rms(src_mono_norm, win_sizes=[win_len])

ifunc = interp1d(linspace(0, len(src_mono_norm)/float(sr), len(w_rms[0])), w_rms[0])
interp_rms = ifunc(times)
centroid_rms = array(centroid[7:]) * interp_rms[7:]
plot(times[7:], centroid_rms, alpha=0.5) #had some 'nan' values in there
title('Spectral Centroid Scaled by RMS')
mean_sc = sum(centroid[7:]) / len(centroid[7:])

print 'Parliament - P-Funk (Wants to Get Funked Up)'
print 'Mean Spectral Centroid: {:.2f} Hz'.format(mean_sc)

src = sources[2][1]
sr = sources[2][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=2048, Fs=sr, window=window_hanning, noverlap=1024)
X = sqrt(Pxx)
centroid = []

for spec in X.T:
    sc = sum(spec*freqs)/sum(spec)
    centroid.append(sc)

win_len = 0.1 * sr
w_rms, win_lens = windowed_rms(src_mono_norm, win_sizes=[win_len])

ifunc = interp1d(linspace(0, len(src_mono_norm)/float(sr), len(w_rms[0])), w_rms[0])
interp_rms = ifunc(times)
centroid_rms = array(centroid) * interp_rms
plot(times, centroid_rms, alpha=0.5)
title('Spectral Centroid Scaled by RMS')
mean_sc = sum(centroid) / len(centroid)

print 'Me - Falling Elevator Blues'
print 'Mean Spectral Centroid: {} Hz'.format(mean_sc)

src = sources[3][1]
sr = sources[3][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=2048, Fs=sr, window=window_hanning, noverlap=1024)
X = sqrt(Pxx)
centroid = []

for spec in X.T:
    sc = sum(spec*freqs)/sum(spec)
    centroid.append(sc)

win_len = 0.1 * sr
w_rms, win_lens = windowed_rms(src_mono_norm, win_sizes=[win_len])

ifunc = interp1d(linspace(0, len(src_mono_norm)/float(sr), len(w_rms[0])), w_rms[0])
interp_rms = ifunc(times)
centroid_rms = array(centroid) * interp_rms
plot(times, centroid_rms, alpha=0.5)
title('Spectral Centroid Scaled by RMS')
mean_sc = sum(centroid) / len(centroid)


print 'James Blake - The Wilhelm Scream'
print 'Mean Spectral Centroid: {} Hz'.format(mean_sc)

src = sources[4][1]
sr = sources[4][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=2048, Fs=sr, window=window_hanning, noverlap=1024)
X = sqrt(Pxx)
centroid = []

for spec in X.T:
    sc = sum(spec*freqs)/sum(spec)
    centroid.append(sc)

win_len = 0.1 * sr
w_rms, win_lens = windowed_rms(src_mono_norm, win_sizes=[win_len])

ifunc = interp1d(linspace(0, len(src_mono_norm)/float(sr), len(w_rms[0])), w_rms[0])
interp_rms = ifunc(times)
centroid_rms = array(centroid[14:]) * interp_rms[14:]# more 'nan' values
plot(times[14:], centroid_rms, alpha=0.5)
title('Spectral Centroid Scaled by RMS')
mean_sc = sum(centroid[14:]) / len(centroid_rms[14:])


print 'Stars of the Lid - Central Texas'
print 'Mean Spectral Centroid: {} Hz'.format(mean_sc)

src = sources[0][1]
sr = sources[0][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=2048, Fs=sr, window=window_hanning, noverlap=1024)
X = sqrt(Pxx)

spread = []

for spec in X.T:
    ss = var(spec)
    spread.append(ss)

print 'Black Keys - Everlasting Light'
plot(spread);
title('Spectral Spread')

src = sources[1][1]
sr = sources[1][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=2048, Fs=sr, window=window_hanning, noverlap=1024)
X = sqrt(Pxx)

spread = []

for spec in X.T:
    ss = var(spec)
    spread.append(ss)

print 'Parliament - P-Funk (Wants to Get Funked Up)'
plot(spread);
title('Spectral Spread')

src = sources[2][1]
sr = sources[2][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=2048, Fs=sr, window=window_hanning, noverlap=1024)
X = sqrt(Pxx)

spread = []

for spec in X.T:
    ss = var(spec)
    spread.append(ss)

print 'Me - Falling Elevator Blues'
plot(spread);
title('Spectral Spread')

src = sources[3][1]
sr = sources[3][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=2048, Fs=sr, window=window_hanning, noverlap=1024)
X = sqrt(Pxx)

spread = []

for spec in X.T:
    ss = var(spec)
    spread.append(ss)

print 'James Blake - The Wilhelm Scream'
plot(spread);
title('Spectral Spread')

src = sources[4][1]
sr = sources[4][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=2048, Fs=sr, window=window_hanning, noverlap=1024)
X = sqrt(Pxx)

spread = []

for spec in X.T:
    ss = var(spec)
    spread.append(ss)

print 'Stars of the Lid - Central Texas'
plot(spread);
title('Spectral Spread')

src = sources[0][1]
sr = sources[0][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=2048, Fs=sr, window=window_hanning, noverlap=1024)
X = sqrt(Pxx)

skewness = []
for spec in X.T:
    ss = moment(spec, moment=3)
    skewness.append(ss)

print 'Black Keys - Everlasting Light'
plot(skewness)
title('Spectral Skewness')

src = sources[1][1]
sr = sources[1][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=2048, Fs=sr, window=window_hanning, noverlap=1024)
X = sqrt(Pxx)

skewness = []
for spec in X.T:
    ss = moment(spec, moment=3)
    skewness.append(ss)

print 'Parliament - P-Funk (Wants to Get Funked Up)'
plot(skewness)
title('Spectral Skewness')

src = sources[2][1]
sr = sources[2][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=2048, Fs=sr, window=window_hanning, noverlap=1024)
X = sqrt(Pxx)

skewness = []
for spec in X.T:
    ss = moment(spec, moment=3)
    skewness.append(ss)

print 'Me - Falling Elevator Blues'
plot(skewness)
title('Spectral Skewness')

src = sources[3][1]
sr = sources[3][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=2048, Fs=sr, window=window_hanning, noverlap=1024)
X = sqrt(Pxx)

skewness = []
for spec in X.T:
    ss = moment(spec, moment=3)
    skewness.append(ss)

print 'James Blake - The Wilhelm Scream'
plot(skewness)
title('Spectral Skewness')

src = sources[4][1]
sr = sources[4][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=4096, Fs=sr, window=window_hanning, noverlap=2048)
X = sqrt(Pxx)

skewness = []
for spec in X.T:
    ss = moment(spec, moment=3)
    skewness.append(ss)

print 'Stars of the Lid - Central Texas'
plot(skewness)
title('Spectral Skewness')

src = sources[0][1]
sr = sources[0][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=4096, Fs=sr, window=window_hanning, noverlap=2048)
X = sqrt(Pxx)

kurtosis = []
for spec in X.T:
    sk = moment(spec, moment=4)
    kurtosis.append(sk)

print 'Black Keys - Everlasting Light'
plot(kurtosis)
title('Spectral Kurtosis')

src = sources[1][1]
sr = sources[1][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=4096, Fs=sr, window=window_hanning, noverlap=2048)
X = sqrt(Pxx)

kurtosis = []
for spec in X.T:
    sk = moment(spec, moment=4)
    kurtosis.append(sk)

print 'Parliament - P-Funk (Wants to Get Funked Up)'
plot(kurtosis)
title('Spectral Kurtosis')

src = sources[2][1]
sr = sources[2][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=4096, Fs=sr, window=window_hanning, noverlap=2048)
X = sqrt(Pxx)

kurtosis = []
for spec in X.T:
    sk = moment(spec, moment=4)
    kurtosis.append(sk)

print 'Me - Falling Elevator Blues'
plot(kurtosis)
title('Spectral Kurtosis')

src = sources[3][1]
sr = sources[3][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=4096, Fs=sr, window=window_hanning, noverlap=2048)
X = sqrt(Pxx)

kurtosis = []
for spec in X.T:
    sk = moment(spec, moment=4)
    kurtosis.append(sk)

print 'James Blake - The Wilhelm Scream'
plot(kurtosis)
title('Spectral Kurtosis')

src = sources[4][1]
sr = sources[4][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=4096, Fs=sr, window=window_hanning, noverlap=2048)
X = sqrt(Pxx)

kurtosis = []
for spec in X.T:
    sk = moment(spec, moment=4)
    kurtosis.append(sk)

print 'Stars of the Lid - Central Texas'
plot(kurtosis)
title('Spectral Kurtosis')

src = sources[0][1]
sr = sources[0][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=1024, Fs=sr, window=window_hanning, noverlap=512)

X_mag = sqrt(Pxx)

X1 = X_mag[:,:-1]/sum(X_mag[:,:-1])
X2 = X_mag[:, 1:]/sum(X_mag[:, 1:])

spec_diff = X2 - X1
flux = sum(spec_diff**2, axis=0)

twinx()

print 'Black Keys - Everlasting Light'
plot(linspace(0, times.max(), len(flux)), flux, lw=1, color='w')
title('Spectral Flux')

src = sources[1][1]
sr = sources[1][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=1024, Fs=sr, window=window_hanning, noverlap=512)

X_mag = sqrt(Pxx)

X1 = X_mag[:,:-1]/sum(X_mag[:,:-1])
X2 = X_mag[:, 1:]/sum(X_mag[:, 1:])

spec_diff = X2 - X1
flux = sum(spec_diff**2, axis=0)

twinx()

print 'Parliament - P-Funk (Wants to Get Funked Up)'
plot(linspace(0, times.max(), len(flux)), flux, lw=1, color='w')
title('Spectral Flux')

src = sources[2][1]
sr = sources[2][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=1024, Fs=sr, window=window_hanning, noverlap=512)

X_mag = sqrt(Pxx)

X1 = X_mag[:,:-1]/sum(X_mag[:,:-1])
X2 = X_mag[:, 1:]/sum(X_mag[:, 1:])

spec_diff = X2 - X1
flux = sum(spec_diff**2, axis=0)

twinx()

print 'Me - Falling Elevator Blues'
plot(linspace(0, times.max(), len(flux)), flux, lw=1, color='w')
title('Spectral Flux')

src = sources[3][1]
sr = sources[3][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=1024, Fs=sr, window=window_hanning, noverlap=512)

X_mag = sqrt(Pxx)

X1 = X_mag[:,:-1]/sum(X_mag[:,:-1])
X2 = X_mag[:, 1:]/sum(X_mag[:, 1:])

spec_diff = X2 - X1
flux = sum(spec_diff**2, axis=0)

twinx()

print 'James Blake - The Wilhelm Scream'
plot(linspace(0, times.max(), len(flux)), flux, lw=1, color='w')
title('Spectral Flux')

src = sources[4][1]
sr = sources[4][0]
# Downmix
src_mono = sum(src.astype(float)/src.shape[1], axis=1).astype(src.dtype)
# Normalize
abs_max = max(abs(src_mono.min().astype(float)), abs(src_mono.max().astype(float)))
src_mono_norm = src_mono.astype(float) / abs_max

Pxx, freqs, times, im = specgram(src_mono_norm, NFFT=1024, Fs=sr, window=window_hanning, noverlap=512)

X_mag = sqrt(Pxx)

X1 = X_mag[:,:-1]/sum(X_mag[:,:-1])
X2 = X_mag[:, 1:]/sum(X_mag[:, 1:])

spec_diff = X2 - X1
flux = sum(spec_diff**2, axis=0)

twinx()

print 'Stars of the Lid - Central Texas'
plot(linspace(0, times.max(), len(flux)), flux, lw=1, color='w')
title('Spectral Flux')