Anonymous avatar Anonymous committed a7a7677

improving stspecgram

Comments (0)

Files changed (4)

 # row ii -> ii/L
 # Then the freq of row ii is
 freq[ii] = ii/Tn
+
+
+------------
+Noter, the NIH imaging center also has multidimensional stockwell transforms
 import stockwell.smt as smt
 # eegfile="/mnt/hgfs/clee/Documents/swaineeg/NKT/EEG2100/CA75510M_1-1+_1-2+.edf"
 eegfile='/home/clee/Dropbox/data/swainAFIB_CA46803E_1-1+.edf'
+eegfile=r'/home/clee/Dropbox/data/swainShortAFIB_CA46803D_1-1+_1-2+.edf'
 open(eegfile)
 e=_edflib.Edfreader(eegfile)
 

stockwell/__init__.py

 """
 from __future__ import division
 
-from st import *
+# import the methods implemented in C
+from st import st, ist, hilbert
 
 def stfreq(f,length, srate):
 	"""

stockwell/plots.py

+"""
+provide matplotlib-based visualization functions for stockwell transforms
+"""
+import stockwell 
 import matplotlib.pyplot as plt
 
 from matplotlib.offsetbox import TextArea, DrawingArea, OffsetImage, \
 
 from mpl_toolkits.axes_grid1 import host_subplot
 
-def plotspec(psx, fs=2.0, lofreq=None, hifreq=None, t1=None, t2=None):
+def plotspec(psx, fs=2.0, lofreq=None, hifreq=None, t0=None, t1=None):
     """
     useful for plotting the power of a stockwell transform
     it relies upon matplotlib for display
     >>> # pylab.show() # to visualize this
     """
     extent = [0,psx.shape[1], 0.0, fs/2.0] # default extents
-    if t1 != None and t2 != None:
-        extent[0] = t1
-        extent[1] = t2
+    if t0 != None and t1 != None:
+        extent[0] = t0
+        extent[1] = t1
     if lofreq != None:
         extent[2] = lofreq
     if hifreq != None:
         extent[3] = hifreq
-    
+    plt.ylabel('Hz')
     return plt.imshow(psx, extent=extent, aspect='auto', origin='lower')
 
 
+def stspecgram(x,fs,lofreq=None, hifreq=None, t0=None, t1=None):
+    """
+    plot out the stockwell spectrum abs(st(x))
+    given frequency sampling fs in Hz
+
+    lofreq and hifreq are the frequency limits expressed in terms of the nyquist frequency(?)
+    """
+
+    n = x.shape[0]
+    if t0==None:
+        t0=0.0
+    if t1==None:
+        t1=n/float(fs)+t0
+
+    if lofreq==None and hifreq==None:
+        sx=stockwell.st(x)
+        return plotspec(abs(sx),fs, t0=t0,t1=t1)
+
+    lorow=stockwell.stfreq(lofreq,n,fs)
+    hirow=stockwell.stfreq(hifreq,n,fs)
+    sx=stockwell.st(x, lorow,hirow)
+    return plotspec(abs(sx), fs, lofreq=lofreq,hifreq=hifreq, t0=t0,t1=t1)
+
 
 def stpowerstack(x,stx):
     """
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.