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.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
import PyGFT as gft
Start with ordinary STFT, so we can compare to something...
In [3]:
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 [4]:
speech = np.loadtxt('benchmark_signals/speech.txt')
In [5]:
tf(speech, 8000)
In [6]:
# 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]
Out[6]:
In [7]:
# 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:
In [8]:
# Get the partitions
partitions = gft.partitions(256)
# Get the windows
windows = gft.gaussianWindows(N)
In [9]:
partitions
Out[9]:
What are these weird numbers?
In [10]:
34359738372/8589934593., 137438953488/34359738372.
Out[10]:
OK, fine. Let's try this...
In [11]:
8589934593/256**4.
Out[11]:
In [12]:
for p in partitions: print p/partitions[-1]
I have no idea what that is...
The windows seem more manageable:
In [13]:
windows.shape
Out[13]:
In [14]:
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()
Finally, interpolate the GFT spectrum and plot a spectrogram
In [15]:
plt.figure()
plt.imshow(abs(gft.gft1dInterpolateNN(SIGNAL)))
plt.show()
Following this document.
In [16]:
syn = np.loadtxt('benchmark_signals/synthetic.txt')
SYN = gft.gft1d(syn, 'gaussian')
SYN
Out[16]:
In [17]:
SYN.shape
Out[17]:
In [18]:
plt.figure(figsize=(12,6))
plt.imshow(abs(gft.gft1dInterpolateNN(SYN)))
plt.show()
In [19]:
tf(syn, 800)
In [20]:
micro = np.loadtxt('benchmark_signals/microseismic.txt')
MICRO = gft.gft1d(micro, 'gaussian')
In [21]:
plt.figure()
plt.imshow(abs(gft.gft1dInterpolateNN(MICRO)))
plt.show()
In [22]:
tf(micro, 2000, w=32)
In [23]:
ls benchmark_signals/
In [24]:
# Crashes Python, something to do with fftw
seismic = np.loadtxt('benchmark_signals/seismic_trace.txt')
#SEISMIC = gft.gft1d(seismic, 'gaussian')
#SEISMIC
In [ ]:
plt.figure()
plt.imshow(abs(gft.gft1dInterpolateNN(SEISMIC)))
plt.show()
In [26]:
seismic[:10]
Out[26]:
In [27]:
tremor = np.loadtxt('benchmark_signals/tremor.txt')
TREMOR = gft.gft1d(tremor, 'gaussian')
TREMOR
Out[27]:
In [28]:
# Crashes fftw?
plt.figure()
#plt.imshow(abs(gft.gft1dInterpolateNN(TREMOR)))
plt.show()
In [29]:
# Crashes Python
tohoku = np.loadtxt('benchmark_signals/tohoku.txt')
#TOHOKU = gft.gft1d(tohoku, 'gaussian')
In [ ]:
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..."
In [38]:
import scipy.signal
In [39]:
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')
In [40]:
ch.shape
Out[40]:
In [42]:
CH = gft.gft1d(ch, 'gaussian')
In [ ]:
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()
In [35]:
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()
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.
In [36]:
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()
In [37]:
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()
In [37]: