Fun with spectrograms

In [1]:

    

import numpy as np
from numpy.fft import fft, fftfreq
import matplotlib.pyplot as plt
%matplotlib inline


    

In [2]:

    

from pylab import specgram

dt = 0.001
t = np.arange(0.0, 3.0, dt)
s1 = np.sin(2*np.pi*20*t)
s2 = 2*np.sin(2*np.pi*20*t)

# create a transient "chirp"
mask = np.where(np.logical_and(t>1, t<2), 1.0, 0.0)
s2 = s2 * mask

# add some noise into the mix
nse = 0.01*np.random.randn(len(t))

x = s1 + s2 + nse # the signal
NFFT = 256      # the length of the windowing segments
Fs = int(1.0/dt)  # the sampling frequency

# Pxx is the segments x freqs array of instantaneous power, freqs is
# the frequency vector, bins are the centers of the time bins in which
# the power is computed, and im is the matplotlib.image.AxesImage
# instance

plt.figure(figsize=(12,8))
plt.subplot(211)
plt.plot(t, x)
ax1 = plt.subplot(212)
Pxx, freqs, bins, im = plt.specgram(x, NFFT=NFFT, Fs=Fs, noverlap=255)
ax1.set_ylim((0,256))
plt.show()


    

    




The history saving thread hit an unexpected error (OperationalError('database is locked',)).History will not be written to the database.






    

In [3]:

    

# Import the libraries we'll need
import urllib2
import StringIO
import zipfile
 
# Read the data from the web.
web_data = urllib2.urlopen('http://www.thesoundarchive.com/simpsons/homer/mcrumble.wav')
 
# Convert the bytes into a file-like object in memory
wav_in_memory = StringIO.StringIO(web_data.read())


    

In [4]:

    

import scipy.io.wavfile as wavfile

rate,data = wavfile.read(wav_in_memory)

plt.figure(figsize=(12,8))
plt.subplot(411)
plt.plot(range(len(data)),data)
plt.subplot(412)

# NFFT is the number of data points used in each block for the FFT
# and noverlap is the number of points of overlap between blocks
plt.specgram(data, NFFT=128, noverlap=64) # small window
plt.subplot(413)
plt.specgram(data, NFFT=512, noverlap=128) 
plt.subplot(414)
plt.specgram(data, NFFT=1024, noverlap=512) # big window

plt.show()


    

    




/Library/Python/2.7/site-packages/scipy-0.14.0.dev_7cefb25-py2.7-macosx-10.9-intel.egg/scipy/io/wavfile.py:172: WavFileWarning: Chunk (non-data) not understood, skipping it.
  WavFileWarning)






    

In [5]:

    

# create a wave with 1Mhz and 0.5Mhz frequencies
dt = 0.001 # 1 ms
t = np.arange(0, 20, dt)
fscale = t/max(t)
y = np.cos(2 * np.pi * 30 * t*fscale) + (np.cos(2 * np.pi * 100 * t*fscale) * np.cos(2 * np.pi * 50 * t*fscale))
y *= np.hanning(len(y))
yy = np.concatenate((y, ([0] * 10 * len(y))))

# FFT of this
Fs = 1 / dt  # sampling rate
n = len(yy)  # length of the signal
k = np.arange(n)
T = n / Fs
frq = k / T  # two sides frequency range
frq = frq[range(n / 2)]  # one side frequency range
Y = fft(yy) / n  # fft computing and normalization
Y = Y[range(n / 2)] / max(Y[range(n / 2)])

plt.figure(figsize=(12,8))

# plotting the data
plt.subplot(3, 1, 1)
plt.plot(t * 1e3, y, 'r')
plt.xlabel('Time (ms)')
plt.ylabel('Amplitude')
plt.grid()

# plotting the spectrum
plt.subplot(3, 1, 2)
plt.plot(frq[0:600], abs(Y[0:600]), 'k')
plt.xlabel('Freq (Hz)')
plt.ylabel('|Y(freq)|')
plt.grid()

# plotting the specgram
ax1 = plt.subplot(3, 1, 3)
Pxx, freqs, bins, im = plt.specgram(y, NFFT=512, Fs=Fs, noverlap=10)
#ax1.set_ylim((0,250))
plt.show()


    

In [6]:

    

import scipy.signal


    

In [7]:

    

timebase = np.arange(0,6,0.001)
f0 = 1.0
t1 = 6.0
f1 = 64.0
ch = scipy.signal.chirp(timebase, f0, t1, f1, method='quadratic')


    

In [8]:

    

plt.plot(ch)
plt.show()


    

In [9]:

    

ax1 = plt.subplot(111)
plt.specgram(ch, Fs=1000)[3]
ax1.set_ylim((0,64))
plt.show()


    

In [10]:

    

f02 = 128.0
t12 = 6.0
f12 = 8.0
ch2 = scipy.signal.chirp(timebase, f02, t12, f12, method='logarithmic')


    

In [11]:

    

ax1 = plt.subplot(111)
plt.specgram(ch+ch2, Fs=1000)[3]
ax1.set_ylim((0,128))
plt.show()


    

In [12]:

    

speech_filename = 'benchmark_signals/speech.txt'
with open(speech_filename) as f:
    s = f.readlines()
s[:10]


    

    
Out[12]:





['  -1.5500000e+02\r\n',
 '  -4.5000000e+01\r\n',
 '   1.3900000e+02\r\n',
 '   1.7100000e+02\r\n',
 '   5.7000000e+01\r\n',
 '  -3.3000000e+01\r\n',
 '  -4.0000000e+00\r\n',
 '   1.5500000e+02\r\n',
 '   2.1100000e+02\r\n',
 '   2.1100000e+02\r\n']

In [13]:

    

speech = np.loadtxt(speech_filename)


    

In [14]:

    

speech.size


    

    
Out[14]:





1400

In [15]:

    

fs_speech = 8000 # sampling frequency, in Hz


    

In [16]:

    

# Let's make a function to plot a signal and its specgram

def tf(signal, fs, w=256, ylim=None):
    '''
    Plot a 1D signal and its spectrogram
    
    signal: a 1D array of amplitudes
    fs: the sampling frequency
    w: the window length
    '''
    
    dt = 1./fs
    n = signal.size
    t = np.arange(0.0, n*dt, dt)
        
    # add some noise into the mix
    #nse = 0.01*np.random.randn(len(t))
    
    noverlap = w - 1
        
    # Pxx is the segments x freqs array of instantaneous power, freqs is
    # the frequency vector, bins are the centers of the time bins in which
    # the power is computed, and im is the matplotlib.image.AxesImage
    # instance
    
    plt.figure(figsize=(12,8))
    plt.subplot(211)
    plt.plot(t, signal)
    ax1 = plt.subplot(212)
    Pxx, freqs, bins, im = plt.specgram(signal, NFFT=w, Fs=fs, noverlap=noverlap, cmap='cubehelix_r')
    if ylim:
        ax1.set_ylim((0,ylim))
    plt.show()


    

In [17]:

    

tf(speech, 8000)


    

In [18]:

    

synthetic = np.loadtxt('benchmark_signals/synthetic.txt')


    

In [19]:

    

tf(synthetic, 800, w=128, ylim=256)


    

    




/Library/Python/2.7/site-packages/matplotlib-override/matplotlib/axes.py:8936: RuntimeWarning: divide by zero encountered in log10
  Z = 10. * np.log10(Pxx)






    

In [20]:

    

seismic = np.loadtxt('benchmark_signals/seismic_trace.txt')


    

In [21]:

    

tf(seismic, 500, w=128, ylim=128)


    

In [22]:

    

tremor = np.loadtxt('benchmark_signals/tremor.txt')


    

In [23]:

    

tf(tremor, 100)


    

In [24]:

    

microseismic = np.loadtxt('benchmark_signals/microseismic.txt')


    

In [25]:

    

tf(microseismic, 2000, w=32)


    

In [26]:

    

tohoku = np.loadtxt('benchmark_signals/tohoku.txt')


    

In [27]:

    

tf(tohoku, 20, w=1024, ylim=0.5)