In [1]:
import numpy as np
%matplotlib nbagg
import matplotlib.pyplot as plt
import scipy.io.wavfile
In [2]:
import os
print os.getcwd()
In [3]:
#https://github.com/chichilalescu/python-for-scientific-computing/tree/master/data
r, d0 = scipy.io.wavfile.read('data/tst0.wav')
print type(d0)
print r
print d0.shape, d0.dtype
t0 = np.array(range(d0.shape[0])).astype(np.float) / r
print t0
r, d1 = scipy.io.wavfile.read('data/tst1.wav')
print r
print d1.shape, d1.dtype
t1 = np.array(range(d1.shape[0])).astype(np.float) / r
print t1
In [4]:
ax = plt.figure().add_subplot(111)
ax.plot(t1, d1[:, 0])
ax.plot(t1, d1[:, 1])
Out[4]:
In [5]:
sfull1 = np.fft.rfft(d1, axis = 0)
sfull0 = np.fft.rfft(d0, axis = 0)
k0 = np.array(range(sfull0.shape[0])).astype(t0.dtype) / t0[-1]
k1 = np.array(range(sfull1.shape[0])).astype(t1.dtype) / t1[-1]
ax = plt.figure(figsize=(6, 4)).add_axes([.15, .15, .7, .7])
ax.plot(k0, np.abs(sfull0)[:, 0]**2, label = '0')
ax.plot(k1, np.abs(sfull1)[:, 0]**2, label = '1')
ax.set_yscale('log')
ax.set_xscale('log')
ax.set_ylabel('intensity')
ax.set_xlabel('Hz')
ax.legend(loc = 'best')
Out[5]:
but we don't want $10^{-1}$ Hz, we only want stuff above 20Hz. or 16, if you have really good ears.
16 Hz means 1/16 seconds, and we have a sampling rate of 44100. Therefore we want 44100/16 samples at most in a signal
In [6]:
print t0.shape[0]/(44100/16.), t1.shape[0]/(44100/16.)
In [7]:
nsamples = int(t0.shape[0]/(44100/16.))
nsamplestmp = int(t1.shape[0]/(44100/16.))
if nsamplestmp < nsamples:
nsamples = nsamplestmp
print nsamples
In [8]:
d0 = d0[:nsamples*int(44100/16.)].reshape(nsamples, int(44100/16.), 2)
d1 = d1[:nsamples*int(44100/16.)].reshape(nsamples, int(44100/16.), 2)
t0 = t0[:int(44100/16.)]
t1 = t1[:int(44100/16.)]
In [9]:
sfull1 = np.fft.rfft(d1, axis = 1)
sfull0 = np.fft.rfft(d0, axis = 1)
k0 = np.array(range(sfull0.shape[1])).astype(t0.dtype) / t0[-1]
print k0
ax = plt.figure(figsize=(6, 4)).add_axes([.15, .15, .7, .7])
s0 = np.average(np.abs(sfull0)[:, :, 0]**2, axis = 0)
s1 = np.average(np.abs(sfull1)[:, :, 0]**2, axis = 0)
ax.plot(k0, s0, label = '0')
ax.plot(k0, s1, label = '1')
ax.set_yscale('log')
ax.set_xscale('log')
ax.set_ylabel('intensity')
ax.set_xlabel('Hz')
ax.legend(loc = 'best')
Out[9]:
In [10]:
ax = plt.figure(figsize=(6, 4)).add_axes([.15, .15, .7, .7])
ax.hist(d0[:, :, 0].ravel(),
bins = 32,
histtype = 'step')
ax.hist(d0[:, :, 1].ravel(),
bins = 32,
histtype = 'step')
ax.hist(d1[:, :, 0].ravel(),
bins = 32,
histtype = 'step')
ax.hist(d1[:, :, 1].ravel(),
bins = 32,
histtype = 'step')
Out[10]:
In [10]: