The fast S-transform, implemented by pyfst
from UofC, as part of PyGFT.
This package implements the fast frequency domain GFT algorithm published by Brown et. al. (2010), in both 1 dimension (signals) and two dimensions (images). The fast GFT is an O(N log N) algorithm, which is the same order as the fast Fourier transform. In general, the GFT implemented in this library will take roughly twice as long to calculate as the FFT of the same signal.
It depends on fftw
, which must be installed first.
Brown, RA, Lauzon, ML, and Frayne, R (2010). A General Description of Linear Time-Frequency Transforms and Formulation of a Fast, Invertible Transform That Samples the Continuous S-Transform Spectrum Nonredundantly. IEEE Trans. Signal Process. 58 (1), pp. 281-290, Jan 2010.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import PyGFT as gft
Start with ordinary STFT, so we can compare to something...
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()
speech = np.loadtxt('benchmark_signals/speech.txt')
tf(speech, 8000)
Here's an attempt at the one in examples.py
:
# Create a signal
N = 256
x = np.arange(0,1,1./N)
signal = np.zeros(len(x), np.complex128)
# Make our signal 1 in the middle, zero elsewhere
signal[N/2] = 1.0
# Do the GFT on the signal
SIGNAL = gft.gft1d(signal, 'gaussian')
SIGNAL[:5]
array([ 0.00012227 +6.50521303e-19j, 0.00012197 -1.20133346e-05j, 0.00012108 -2.40845559e-05j, 0.00011957 -3.62726693e-05j, 0.00011742 -4.86389636e-05j])
# Plot the magnitude of the signal
plt.figure()
ax1 = plt.subplot(211)
ax1.plot(x, signal.real)
ax2 = plt.subplot(212)
ax2.plot(x, SIGNAL.real)
ax2.set_xlim((0,1))
plt.show()
Visualize the partitions and the windows:
# Get the partitions
partitions = gft.partitions(256)
# Get the windows
windows = gft.gaussianWindows(N)
partitions
array([ 8589934593, 34359738372, 137438953488, 549755813952, 962072674496, 1065151889648, 1090921693436, 1099511628031, 4294967295])
What are these weird numbers?
34359738372/8589934593., 137438953488/34359738372.
(4.0, 4.0)
OK, fine. Let's try this...
8589934593/256**4.
2.0000000002328306
for p in partitions: print p/partitions[-1]
2 8 32 128 224 248 254 256 1
I have no idea what that is...
The windows seem more manageable:
windows.shape
(256,)
plt.figure(figsize=(15, 6))
ax1 = plt.subplot(111)
ax1.plot(x, SIGNAL.real)
# For each partition, draw a line on the graph. Since the partitions only indicate the
# positive frequencies, we need to draw partitions on both sides of the DC
for i in partitions:
print 0.5 + (i/N**4.)/float(N)
print 0.5 - (i/N**4.)/float(N)
ax1.axvline(0.5+(i/N**4.)/float(N), color='r', alpha=0.2)
ax1.axvline(0.5-(i/N**4.)/float(N), color='r', alpha=0.2)
# Plot the windows
w = gft.shift(windows, 128)
ax1.plot(x, w)
ax1.set_xlim((0,1))
plt.show()
0.507812500001 0.492187499999 0.531250000004 0.468749999996 0.625000000015 0.374999999985 1.00000000006 -5.82076609135e-11 1.37500000017 -0.375000000175 1.46875000022 -0.468750000218 1.49218750023 -0.492187500229 1.50000000023 -0.500000000232 0.503906249999 0.496093750001
Finally, interpolate the GFT spectrum and plot a spectrogram
plt.figure()
plt.imshow(abs(gft.gft1dInterpolateNN(SIGNAL)))
plt.show()
/usr/lib/pymodules/python2.7/matplotlib/colors.py:533: RuntimeWarning: invalid value encountered in less cbook._putmask(xa, xa<0.0, -1)
Following this document.
syn = np.loadtxt('benchmark_signals/synthetic.txt')
SYN = gft.gft1d(syn, 'gaussian')
SYN
array([ 0.40429712+0.17552105j, -0.09675209-0.05918993j, 0.00072977-0.00167102j, ..., -0.00688487+0.00143456j, -0.00353827+0.00164455j, -0.19792998+0.05953251j])
SYN.shape
(8192,)
plt.figure(figsize=(12,6))
plt.imshow(abs(gft.gft1dInterpolateNN(SYN)))
plt.show()
tf(syn, 800)
micro = np.loadtxt('benchmark_signals/microseismic.txt')
MICRO = gft.gft1d(micro, 'gaussian')
plt.figure()
plt.imshow(abs(gft.gft1dInterpolateNN(MICRO)))
plt.show()
tf(micro, 2000, w=32)
ls benchmark_signals/
microseismic.txt* speech.txt* tohoku.txt* seismic_trace.txt* synthetic.txt* tremor.txt*
# Crashes Python, something to do with fftw
seismic = np.loadtxt('benchmark_signals/seismic_trace.txt')
#SEISMIC = gft.gft1d(seismic, 'gaussian')
#SEISMIC
plt.figure()
plt.imshow(abs(gft.gft1dInterpolateNN(SEISMIC)))
plt.show()
seismic[:10]
array([ 6754.2148, 10132.801 , 7748.4023, 4086.103 , 2672.2512, 1717.0735, -2029.7522, -7296.2305, -8909.8672, -4533.5469])
tremor = np.loadtxt('benchmark_signals/tremor.txt')
TREMOR = gft.gft1d(tremor, 'gaussian')
TREMOR
array([ 102.43758809 -3.35443605j, -68.65726164 +8.91291227j, -9.93064552 +49.85290072j, ..., -111.46516574 +98.54789685j, 112.14116362-135.39792285j, -67.10028747 +88.20753371j])
# Crashes fftw?
plt.figure()
#plt.imshow(abs(gft.gft1dInterpolateNN(TREMOR)))
plt.show()
<matplotlib.figure.Figure at 0x7f55fc51f3d0>
# Crashes Python
tohoku = np.loadtxt('benchmark_signals/tohoku.txt')
#TOHOKU = gft.gft1d(tohoku, 'gaussian')
plt.figure()
plt.imshow(abs(gft.gft1dInterpolateNN(TOHOKU)))
plt.show()
For whatever reason, this crashes IPython on one of my machines. "The kernel appears to have died..."
import scipy.signal
timebase = np.arange(0,8.192,0.004)
f0 = 1.0
t1 = 6.0
f1 = 64.0
ch = scipy.signal.chirp(timebase, f0, t1, f1, method='quadratic')
ch.shape
(2048,)
CH = gft.gft1d(ch, 'gaussian')
plt.figure(figsize=(15,8))
ax1 = plt.subplot(111)
# Sometimes crashes
interpol = abs(gft.gft1dInterpolateNN(CH))[:ch.size/2]
interpol = np.flipud(interpol)
ax1.pcolor(20*np.log10(interpol), cmap="Greys")
plt.show()
timebase = np.arange(0, 1.024, 0.004)
f1 = 10
f2 = 40
s = 0.2*np.sin(2*np.pi*f1*timebase) + 0.1*np.sin(2*np.pi*f2*timebase)
S = gft.gft1d(s)
plt.figure(figsize=(10,10))
ax1 = plt.subplot(211)
ax1.plot(timebase, s)
ax2 = plt.subplot(212)
interpol = abs(gft.gft1dInterpolateNN(S))[:s.size/2]
interpol = np.flipud(interpol)
ax2.imshow(20*np.log10(interpol), cmap="Greys")
plt.show()
-c:13: RuntimeWarning: divide by zero encountered in log10
Let's try to reproduce some of the figures from Naghizadeh & Innanen (2011).
Mostafa Naghizadeh and Kris A. Innanen (2011). Interpolation of nonstationary seismic records using a fast non-redundant S-transform. SEG Annual Meeting, San Antonio. Available online.
n = 256
#t = np.arange(0,1,1./n)
t = np.arange(0,n)
s = np.zeros(len(t), np.complex128)
s[111:142] = 1.
plt.figure(figsize=(15,3))
ax1 = plt.subplot(111)
ax1.plot(t, s.real, '.')
ax1.plot(t, s.real, lw=0.5)
ax1.set_xlim((0, 256))
ax1.set_ylim((-0.1, 1.1))
plt.xlabel('Time samples')
plt.ylabel('Amplitude')
plt.show()
S = gft.gft1d(s)
tfr = abs(gft.gft1dInterpolateNN(S))
plt.figure(figsize=(12,8))
plt.imshow(tfr[:128,], cmap='Greys', vmax=10)
plt.colorbar(shrink=0.5)
plt.show()