In [1]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

import numpy as np
from librosa.core import stft, istft

from play import play


    

    




/usr/local/lib/python2.7/dist-packages/librosa/core/audio.py:37: UserWarning: Could not import scikits.samplerate. Falling back to scipy.signal
  warnings.warn('Could not import scikits.samplerate. '

In [15]:

    

# freq = 1981
freq = 2000
gensin = np.sin(np.linspace(0, freq*np.pi, 44100))
# gensin *= (2**15-1)


    

In [16]:

    

# gensin = np.array(gensin, dtype='int16')


    

In [17]:

    

plot(gensin[:2000])


    

    
Out[17]:





[<matplotlib.lines.Line2D at 0x7f07e60caa50>]






    

In [6]:

    

play(gensin)


    

In [18]:

    

dsin = stft(gensin, n_fft=2048, hop_length=2048)


    

In [19]:

    

def imshowsq(m):
    """ A helper for showing spectrograms forces to square aspect ratio with no interpolation """
    imshow(m, aspect=float(m.shape[1]) / m.shape[0], interpolation='none')
    colorbar()


    

In [20]:

    

binlo = 38
binhi = 54


    

In [21]:

    

imshowsq(dsin.real[binlo:binhi,:])
title('bins {} through {} of real part of frequency domain'.format(binlo, binhi))


    

    
Out[21]:





<matplotlib.text.Text at 0x7f07e6098690>






    

In [22]:

    

imshowsq(dsin.imag[binlo:binhi,:])
title('bins {} through {} of imaginary part of frequency domain'.format(binlo, binhi))


    

    
Out[22]:





<matplotlib.text.Text at 0x7f07e5efd110>






    

In [108]:

    

plot(dsin.imag[binlo:binhi,10])
plot(dsin.imag[binlo:binhi,11])
plot(dsin.imag[binlo:binhi,12])
plot(dsin.imag[binlo:binhi,13])
title('a few slices of imaginary space at a few particular times')


    

    
Out[108]:





<matplotlib.text.Text at 0x7fbc697ebe50>






    

In [113]:

    

imshowsq(abs(dsin)[binlo:binhi,:])
title('absolute value (maginude complex vector) of frequency domain')


    

    
Out[113]:





<matplotlib.text.Text at 0x7fbc6933ae10>






    

In [114]:

    

imshowsq(np.angle(dsin)[binlo:binhi,:])
title('angle of frequency domain')


    

    
Out[114]:





<matplotlib.text.Text at 0x7fbc6920ab10>






    

In [118]:

    

plot(np.abs(dsin)[38:54,11])
title('slice in time of magnitude')


    

    
Out[118]:





<matplotlib.text.Text at 0x7fbc68e505d0>






    

In [119]:

    

plot(np.angle(dsin)[38:54,11])
title('same time slice of angle')


    

    
Out[119]:





<matplotlib.text.Text at 0x7fbc68d72f50>






    

In [126]:

    

plot(np.sum(abs(dsin), axis=1))
title('sum of magnitudes over time')
xlim(binlo-20,binhi+20)


    

    
Out[126]:





(18, 74)






    

In [128]:

    

plot(np.sum(dsin.real, axis=1))
title('sum of real over time')
xlim(binlo-20,binhi+20)


    

    
Out[128]:





(18, 74)






    

In [131]:

    

# this appears to be some beat frequency combining sin frequency and fft window frequency
plot(dsin.imag[46,:])
title('imaginary value at target frequency bin over time')


    

    
Out[131]:





<matplotlib.text.Text at 0x7fbc68630c90>






    

In [132]:

    

plot(dsin.real[46,:])
title('real value at target frequency bin over time')


    

    
Out[132]:





<matplotlib.text.Text at 0x7fbc683ec910>






    

In [133]:

    

plot(np.abs(dsin[46,:]))
title('maginude at target frequency bin over time')


    

    
Out[133]:





<matplotlib.text.Text at 0x7fbc6837fe10>






    

In [134]:

    

plot(np.angle(dsin)[46,:])
title('angle at target frequency bin over time')


    

    
Out[134]:





<matplotlib.text.Text at 0x7fbc6855cc10>






    

In [78]:

    

np.fft.fftfreq(dsin.shape[0])*44100


    

    
Out[78]:





array([   0.        ,   43.02439024,   86.04878049, ..., -129.07317073,
        -86.04878049,  -43.02439024])

In [79]:

    

# given number of bins in freq space in D and sample rate, lookup
# frequence of sound in bin 46 which matches ~2khz sin wave
# generated above
(np.fft.fftfreq(dsin.shape[0])*44100)[46]


    

    
Out[79]:





1979.1219512195121

In [30]:

    

# TODO: now show same plots, but with a different frequency band closer to 1979
# TODO: any way to make these plots change as a slider changes the frequency