%pylab inline
rcParams['figure.figsize'] = (10, 4) #wide graphs by default
from __future__ import print_function
from __future__ import division
Populating the interactive namespace from numpy and matplotlib
http://en.wikipedia.org/wiki/Variance
For equally likely values:
$$\operatorname{Var}(X) = \frac{1}{n} \sum_{i=1}^n (x_i - \mu)^2$$samps = normal(scale=2.0, size= 10000)
samps.var()
4.0424222584326799
samps.std(), samps.std()**2
(2.0105775932384904, 4.0424222584326808)
a = rand(100)
b = rand(100)
cov(a,b)
array([[ 0.07082438, -0.00174744], [-0.00174744, 0.08639426]])
a = rand(100)* 100
b = rand(100) * 100
cov(a,b)
array([[ 717.56804051, 39.28797053], [ 39.28797053, 944.18916769]])
a = rand(100000)* 100
b = rand(100000) * 100
cov(a,b)
array([[ 832.27860929, -1.08306774], [ -1.08306774, 833.92500888]])
var(a), var(b)
(832.27028649959755, 833.91666963464957)
a = rand(100)
b = rand(100)
c = rand(100)
cov(a + (2*c),b + (2*c))
array([[ 0.35149636, 0.28129156], [ 0.28129156, 0.3858564 ]])
a = rand(100000)
b = rand(100000)
c = rand(100000)
cov(a + (2*c),b + (2*c))
array([[ 0.41751468, 0.33404553], [ 0.33404553, 0.41662556]])
a = rand(100)
b = rand(100)
d = sin(linspace(0, 2*pi*7, 100))
plot(a+d)
plot(b+d)
[<matplotlib.lines.Line2D at 0x7fcd848cbb90>]
cov(a + d,b + d)
array([[ 0.53132808, 0.46277415], [ 0.46277415, 0.57308237]])
plot(a+d + 2)
plot(-2*(b+d))
cov(a + d + 2,-2*(b+d))
array([[ 0.53132808, -0.92554829], [-0.92554829, 2.29232948]])
plot(a + (d*3))
plot(-(b + (d*3)))
[<matplotlib.lines.Line2D at 0x7fcd82cea2d0>]
cov(a + (d*3),-(b + (d*3)))
array([[ 4.43271868, -4.40869777], [-4.40869777, 4.56353902]])
sig1 = sin(linspace(0, 100*2*pi, 44100))
sig2 = sin(linspace(0.1, 0.1 + (100*2*pi), 44100))
plot(sig1)
plot(sig2)
xlim((0, 1000))
(0, 1000)
lags, c, lines, line = xcorr(sig1, sig2, maxlags=500, usevlines=False);
grid()
lags, c, lines, line = xcorr(sig1, sig2, maxlags=50, usevlines=False);
grid()
len(c)
101
plot(c)
[<matplotlib.lines.Line2D at 0x7fcd829d6190>]
argmax(c)
57
cc_index = lags[argmax(c)]
print(cc_index)
7
plot(lags, c)
vlines(cc_index, 0.5, 1.09, color='r', lw=2)
text(cc_index + 1, 1.05, 'CC peak', color='r')
grid()
sr = 44100
f = 100.0
Ps = sr/f
Ps # period in number of samples
441.0
cc_index
7
CC can be used to determine pitch and periodicity.
cc_index/Ps
0.015873015873015872
2 * pi * cc_index/Ps
0.09973310011396169
2 * pi * (cc_index + 1)/Ps
0.11398068584452764
Two different noise samples produce low values of cross-correlation:
noise1 = random.random(44100) - 0.5
noise2 = random.random(44100) - 0.5
lags, c, lines, line = xcorr(noise1, noise2, maxlags=500, usevlines=False);
grid()
And no clear patterns in the CC function:
noise1 = normal(size=44100)
noise2 = normal(size=44100)
lags, c, lines, line = xcorr(noise1, noise2, maxlags=500, usevlines=False);
grid()
Cross-correlation can find harder to see relationships:
a = (rand(44100)*2) - 1
plot(a)
[<matplotlib.lines.Line2D at 0x7fcd82851c50>]
b = sin(linspace(0, 2*pi*200, 44100))*0.3
plot(a+b)
[<matplotlib.lines.Line2D at 0x7fcd82c08910>]
c = (rand(44100)*2) - 1
d = sin(linspace(0.7, 0.7 + (2*pi*200), 44100))*0.3
plot(c+d)
[<matplotlib.lines.Line2D at 0x7fcd82da5a90>]
plot(c+d)
xlim(0,1000)
(0, 1000)
lags, c, lines, line = xcorr(a+b, c+d, maxlags=500);
m = argmax(c[500:600]) + 500
lags[m]
21
sr = 44100
f = 200.0
Ps = sr/f
Ps # period in number of samples
220.5
lags[m]/Ps
0.095238095238095233
(lags[m]/Ps) * 2 * pi
0.59839860068377004
lags, c, lines, line = acorr(a + b, maxlags=500);
plot(c[501:])
vlines(220.5, -0.15, 0.15)
<matplotlib.collections.LineCollection at 0x7fcd80d28150>
CC for Inter-aural time delay (or stereo relationships) using running cross-correlation
from scipy.io import wavfile
sr, in_sig = wavfile.read('cars.wav')
win_size = 1024
hop = 256
cc_lags = 0.001*sr # 1 ms max lag
rccf = []
win_start = arange(0, in_sig.shape[0] - win_size, hop)
for start in win_start:
win = in_sig[start: start+win_size].astype(float)
lags, c, lines, line = xcorr(win[:,0], win[:,1], maxlags=cc_lags)
rccf.append(c)
imshow(array(rccf).T, aspect='auto')
<matplotlib.image.AxesImage at 0x7f6791f03d10>
plot(in_sig)
[<matplotlib.lines.Line2D at 0x7f7cbe28ae10>, <matplotlib.lines.Line2D at 0x7f7cbe295150>]
CC for tempo and rhythm estimation:
sr, in_sig = wavfile.read('superstition.wav')
!aplay superstition.wav
Playing WAVE 'superstition.wav' : Signed 16 bit Little Endian, Rate 44100 Hz, Stereo
plot(in_sig)
[<matplotlib.lines.Line2D at 0x7f7cbe159e50>, <matplotlib.lines.Line2D at 0x7f7cbe0e4190>]
in_sig.shape
(437078, 2)
lags,c, lines, line = acorr(in_sig[:100000,0].astype(double), maxlags= 50000);
i = imread('galileo.png')
imshow(i)
<matplotlib.image.AxesImage at 0x7f7cb4dfea10>
plot(i[:,:,3].T);
i.shape
(1213, 780, 4)
i.dtype
dtype('float32')
i = sum(i[:,:,:-1], axis=2)/3.0
imshow(i)
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7f7ce31df5a8>
i.shape
(1213, 780)
imshow(i)
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7f7ce2de1710>
imshow(i, cmap=cm.gray)
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7f7ce317e710>
i = where(i > 0.9, 0, 1)
imshow(i, cmap='gray')
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7f7cd12245a8>
o = imread('o.png')
imshow(o)
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7f7cd111ca28>
o = o.astype(float).sum(axis=-1)/3
o = where(o > 0.9, 0, 1)
imshow(o, cmap=cm.gray, interpolation='nearest')
<matplotlib.image.AxesImage at 0x7f7cd100bc50>
from scipy.signal import correlate2d
cc = correlate2d(i, o)
imshow(cc)
colorbar()
gcf().set_figheight(8)
imshow(cc, cmap=cm.gray)
colorbar()
gcf().set_figheight(8)
imshow(where(cc > 200, 1, 0), interpolation='nearest', cmap=cm.gray)
colorbar()
gcf().set_figheight(8)
imshow(where(cc > 300, 1, 0), interpolation='nearest', cmap=cm.gray)
colorbar()
gcf().set_figheight(8)
from scipy.ndimage.filters import maximum_filter
imshow(maximum_filter(cc, (10,10)))
gcf().set_figheight(8)
subplot(121)
imshow(maximum_filter(cc, (50,50)))
subplot(122)
imshow(i, cmap=cm.gray)
gcf().set_figheight(8)
mf = maximum_filter(cc, (50,50))
argmax(mf)
215823
unravel_index(argmax(mf), mf.shape)
(262, 721)
argmax(cc)
235551
unravel_index(argmax(cc), cc.shape)
(286, 745)
By: Andrés Cabrera mantaraya36@gmail.com For MAT course MAT 201A at UCSB
This ipython notebook is licensed under the CC-BY-NC-SA license: http://creativecommons.org/licenses/by-nc-sa/4.0/