Exploring reassignment with libtfr, a nicely documented library from Dan Meliza, available on GitHub.
In [1]:
import numpy as np
from numpy.fft import fft, fftfreq
import matplotlib.pyplot as plt
%matplotlib inline
from pylab import specgram
Using the Homer simpson voice from the Glowing Python blog.
In [2]:
# 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 [3]:
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()
In [4]:
data = np.array(data, dtype=float)
data
Out[4]:
In [5]:
data.shape
Out[5]:
In [6]:
import libtfr
In [7]:
N = 256
step = 10
NW = 3.5
In [8]:
w = np.hamming(N)
S_s = libtfr.stft(data, w, step)
S_s.shape
Out[8]:
In [9]:
S_up = np.flipud(S_s)
plt.figure(figsize=(18,5))
plt.imshow(20*np.log10(S_up), aspect=10)
plt.colorbar(shrink=0.9)
#ax1.set_ylim((0,64))
plt.show()
Using discrete prolate spherical sequence (DPSS) tapers.
We can also do the 1D transform:
In [10]:
J = libtfr.mtfft(data, NW)
J.shape
Out[10]:
I don't know what this is.
Let's try the spectrogram...
In [11]:
S_m = libtfr.mtm_spec(data, N, step, NW)
S_m.shape
Out[11]:
In [12]:
S_mup = np.flipud(S_m)
plt.figure(figsize=(18,5))
plt.imshow(20*np.log10(S_mup), aspect=10)
plt.colorbar(shrink=0.9)
#ax1.set_ylim((0,64))
plt.show()
In [13]:
# Input signal is real
N = 1024 # number of frequency points
Np = 512 # window size, should be <= N
step = 40 # step size, in time points
K = 6 # number of tapers, default 6
tm = 1.0 # time support of tapers, default 6.0
flock = 0.01 # freq locking; power is not reassigned more than this value, default 0.01
tlock = 5 # time lock, in frames, default 5
# Output freq bins, monotonically increasing,
# 1.0 is Nyquist. Default linear scale with N points.
fgrid = np.linspace(0.1, 0.475, 512)
In [14]:
#S = libtfr.tfr_spec(data, N, shift, Np, K, tm, flock, tlock)
S_t = libtfr.tfr_spec(data, N, step, Np)
S_t.shape
Out[14]:
In [15]:
S_tup = np.flipud(S_t)
plt.figure(figsize=(18,5))
plt.imshow(20*np.log10(S_tup), aspect=0.67)
plt.colorbar(shrink=0.6)
plt.show()
In [16]:
import scipy.signal
In [17]:
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')
f02 = 128.0
t12 = 6.0
f12 = 8.0
ch2 = scipy.signal.chirp(timebase, f02, t12, f12, method='logarithmic')
In [18]:
ax1 = plt.subplot(111)
plt.specgram(ch+ch2, Fs=1000)[3]
ax1.set_ylim((0,128))
plt.show()
In [19]:
N = 2048 # number of frequency points
Np = 512 # window size, should be <= N
step = 40 # step size, in time points
K = 6 # number of tapers, default 6
tm = 6.0 # time support of tapers, default 6.0
chrp = ch + ch2
CHRP_t = libtfr.tfr_spec(chrp, N, step, Np, K, tm)
CHRP_tup = np.flipud(CHRP_t)[-256:,:]
plt.figure(figsize=(8,5))
plt.imshow(20*np.log10(CHRP_tup), aspect=0.4)
plt.colorbar(shrink=0.6)
plt.show()
In [20]:
speech_filename = 'benchmark_signals/speech.txt'
with open(speech_filename) as f:
s = f.readlines()
s[:10]
Out[20]:
In [21]:
speech = np.loadtxt(speech_filename)
In [22]:
speech.size
Out[22]:
In [23]:
fs_speech = 8000 # sampling frequency, in Hz
In [153]:
# Let's make a function to plot a signal and its specgram
def tf(signal, fs, w=256, wtime=False, poverlap=None, ylim=None, colorbar=False, vmin=None, vmax=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))
if wtime:
# Then the window length is time so change to samples
w *= fs
w = int(w)
if poverlap:
# Then overlap is a percentage
noverlap = int(w * poverlap/100.)
else:
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='Greys', vmin=vmin, vmax=vmax)
if colorbar: plt.colorbar()
if ylim:
ax1.set_ylim((0,ylim))
plt.show()
In [147]:
tf(speech, 8000)
In [142]:
synthetic = np.loadtxt('benchmark_signals/synthetic.txt')
In [143]:
tf(synthetic, 800, w=128, ylim=256, vmin=-30)
In [144]:
# Input signal is real
N = 2048 # number of frequency points
Np = 512 # window size, should be <= N
step = 40 # step size, in time points
K = 5 # number of tapers, default 6
tm = 6.0 # time support of tapers, default 6.0
flock = 0.01 # freq locking; power is not reassigned more than this value, default 0.01
tlock = 5 # time lock, in frames, default 5
# Output freq bins, monotonically increasing,
# 1.0 is Nyquist. Default linear scale with N points.
fgrid = np.linspace(0.1, 0.475, 512)
In [59]:
#S = libtfr.tfr_spec(data, N, shift, Np, K, tm, flock, tlock)
SYN_t = libtfr.tfr_spec(synthetic, N, step, Np, K, tm)
SYN_t.shape
Out[59]:
In [97]:
SYN_tup = np.flipud(SYN_t)
plt.figure(figsize=(18,5))
plt.imshow(20*np.log10(SYN_tup), cmap="Greys", aspect=0.1, vmin=-10)
plt.colorbar(shrink=0.6)
plt.show()
In [98]:
seismic = np.loadtxt('benchmark_signals/seismic_trace.txt')
In [103]:
tf(seismic, 500, w=128, ylim=128, vmin=40)
In [33]:
# Input signal is real
N = 1024 # number of frequency points
Np = 128 # window size, should be <= N
step = 20 # step size, in time points
K = 6 # number of tapers, default 6
tm = 6.0 # time support of tapers, default 6.0
flock = 0.01 # freq locking; power is not reassigned more than this value, default 0.01
tlock = 5 # time lock, in frames, default 5
# Output freq bins, monotonically increasing,
# 1.0 is Nyquist. Default linear scale with N points.
fgrid = np.linspace(0.1, 0.475, 512)
In [34]:
#S = libtfr.tfr_spec(data, N, shift, Np, K, tm, flock, tlock)
SEIS_t = libtfr.tfr_spec(seismic, N, step, Np, K, tm)
SEIS_t.shape
Out[34]:
In [108]:
SEIS_tup = np.flipud(SEIS_t)
plt.figure(figsize=(18,5))
plt.imshow(20*np.log10(SEIS_tup), aspect=0.02, cmap="Greys", vmin=140)
plt.colorbar(shrink=0.6)
plt.show()
In [109]:
tremor = np.loadtxt('benchmark_signals/tremor.txt')
In [125]:
tf(tremor, 100, vmin=10)
In [124]:
# Input signal is real
N = 2048 # number of frequency points
Np = 512 # window size, should be <= N
step = 40 # step size, in time points
K = 6 # number of tapers, default 6
tm = 6.0 # time support of tapers, default 6.0
TREMOR_t = libtfr.tfr_spec(tremor, N, step, Np, K, tm)
TREMOR_tup = np.flipud(TREMOR_t)
plt.figure(figsize=(18,5))
plt.imshow(20*np.log10(TREMOR_tup), aspect=0.5, cmap="Greys", vmin=50)
plt.colorbar(shrink=0.6)
plt.show()
In [126]:
microseismic = np.loadtxt('benchmark_signals/microseismic.txt')
In [152]:
tf(microseismic, 2000, w=0.03, wtime=True, poverlap=99, vmin=0)
In [41]:
tohoku = np.loadtxt('benchmark_signals/tohoku.txt')
In [163]:
tf(tohoku, 20, w=100, wtime=True, ylim=0.5, vmin=70)
In [ ]: