stft.pkg

LastWave

#;; -*- Mode: Tcl -*-
#..........................................................................
#                                                                         
#      L a s t W a v e    P a c k a g e 'stft' 2.0
#
#      Author Remi Gribonval
#      
#      File associated to the stft package
#                                                                    
#..........................................................................

# Set the tabulation directory for 
# the short time fourier transforms 
# (and the matching pursuit)
StftTabulationDirectory=_scriptDirectory[0]+'/stft/tabulation'


#
# Dealing with STFT graphs
#

# The default position of stft windows
disp.stft.rect={370 55 600 400}

# Perform zoom on stfts with the mouse
SetZoomBindings GraphStft {'rect' 'xrect'}

#
# Perform cursor displaying on stfts
#

# 
# Regular cursor 
#

setproc _CursorTextGraphStft {}  {
  stft=[setg @object -graph]

  # The 'x' coordinate in sample-coordinate 
  index = int((@x-stft.x0)/stft.dx+.5)
  if (index < 0) {index=0} 
  if (index > (stft.signalSize-1)) {index=stft.signalSize-1} 

  # The 'y' coordinate in sample coordinate
  freq = int(@y*2*stft.freqIdNyquist*stft.dx+.5)
  if (freq < 0) {freq=0} 
  if (freq > stft.freqIdNyquist) {freq=stft.freqIdNyquist} 

  # The time-frequency grid
  {tr tl fr fl} = stft.grid
  T = tr*round(index/tr)
  F = fr*round(freq/fr)
  if ((T>=0) && (T<stft.signalSize) && (F>=0) && (F<=stft.freqIdNyquist) ) {
      a=[stftmax stft -T T T -F F F][0]
  } else {
      a=0
  }
#  stftbest stft atom -T index index -F freq freq]
#  c = [bestchirp atom stft]

  return "$@objname : ($@x,$@y) [id = ($index,$freq)] val = $a"
}

setproc _DrawCursorNoneGraphStft {{&var cursor}} {
  cursor.erase=null
  return [_CursorTextGraphStft]
}

setproc _DrawCursorGraphStft {{&var cursor}}  {
  _ViewDrawCrossHair cursor.view $@x $@y
  cursor.erase=%%`_ViewDrawCrossHair '$cursor.view' $@x $@y`
  return [_CursorTextGraphStft]
}

SetCursorBindings GraphStft {%_DrawCursorNoneGraphStft %_DrawCursorGraphStft}


####################################################
#
# Vertical and Horizontal sections of Stft displays
#
####################################################
  
binding delete 'stftsection'

setproc _StftSection {flagVertical} {
   
  # Get the stft
  stft=[setg $@objname -graph]

  if (flagVertical) {
      # Get the right time index
      index=int((@x-stft.x0)/stft.dx+.5)
      {tr tl fr fl} = stft.grid
      index=tr*int(index/tr+0.5)
      if (index < 0) {index=0}
      if (index >= tr*tl) {index=tr*(tl-1)}
  } else {
      # Get the right freq index
      index=int(@y*stft.dx*2*stft.freqIdNyquist+.5)
      {tr tl fr fl} = stft.grid
      index=fr*int(index/fr+0.5)
      if (index < 0) {index=0}
      if (index >= fr*fl) {index=fr*(fl-1)}
  }
  # If the stft is zero don't display anything
  {max maxt maxf}=[stftmax stft]
  dynamic=[setg $@objname -expo]
  if ($max==0) {
    errorf "_StftSection : the displayed stft is uniformly zero"
  }

  type=stft.type
  if (type=="complex") {
      if (flagVertical) {
	  {coeffR coeffI} = ${index}stft
      } else {
	  {coeffR coeffI} = .${index}stft
      }
      section=coeffR*coeffR+coeffI*coeffI
  } elseif (type=="real") {
      if (flagVertical) {
	  section = ${index}stft
      } else {
	  section = .${index}stft
      }
  } elseif (type=="phase") {
      if (flagVertical) {
	  section = ${index}stft
      } else {
	  section = .${index}stft
      }
# We should unwrap the phase !
  } else {
      errorf "Unnown type %s" type
  }

  if (flagVertical) {
      xLabel="Hz"
  } else {
      xLabel = "seconds"
  }
  # In any case except for a phase display, we must
  # convert to decibel display, relative to the max
  # with the right dynamic range
  if (type!="phase") {
    threshold=((abs(section)/$max)>=10^(-dynamic/10))
    section= threshold*abs(section)/$max+10^(-dynamic/10)*(1-threshold) 
    section= 10*log(section)/log(10) 
    # Then just display the energy
    yLabel="dB"
  } else {
    yLabel="phi"
  }

  # Case the window already exists  
  if ([msge stftWindow exist]) {
      import args 1 &array bindings.GraphStft.StftSection
      updatexaxis=0
      updateyaxis=0
      
      if (StftSection.type!=flagVertical) {
	  updatexaxis=1
	  StftSection.type=flagVertical
      } else {
	  {x0 x1 y0 y1}=[setg stftWindow.fv1.view -bound]
	  {imin imax}=[stats minmax section]
	  ymin=section.Y[imin]
	  ymax=section.Y[imax]
	  if (y0 > ymin) {
	      y0=ymin
	      updateyaxis=1
	  }
	  if (y1 < ymax) {
	      y1=ymax
	      updateyaxis=1
	  }
      }
      if (updatexaxis) {
	  disp stftWindow section
	  setgu ..fv1 -xLabel xLabel
	  setgu ..fv1 -yLabel yLabel
      } elseif (updateyaxis) {
	  setgu stftWindow.fv1 -graph1 section "1" -xLabel xLabel -yLabel yLabel -bound '*' '*' y0 y1
      } else {
	  setgu stftWindow.fv1 -graph1 section "1" -xLabel xLabel -yLabel yLabel
      }
  } else {
      import args 1 &array bindings.GraphStft.StftSection
      StftSection.type=flagVertical
      disp stftWindow section
      setgu ..fv1 -xLabel xLabel
      setgu ..fv1 -yLabel yLabel
  }  
}

setbinding 'stftsection' "{{Shift + Middle button = vertical (decibel) section}} \
{{Ctrl + Middle button = horizontal (decibel) section}}"
setbinding 'stftsection' GraphStft middleButtonDown shift {_StftSection 0}
setbinding 'stftsection' GraphStft middleButtonMotion shift {_StftSection 0}
setbinding 'stftsection' GraphStft middleButtonDown ctrl {_StftSection 1}
setbinding 'stftsection' GraphStft middleButtonMotion ctrl {_StftSection 1}

binding activate 'stftsection'

###############################################
#
# Vertical sections
#
###############################################
  
binding delete 'fftsection'

setproc _FFTSection {} {
   
  # Get the signal
  sig=[setg $@objname -graph]

  # Get the FFT parameters
  import args 1 bindings.GraphSignal.FftSection.windowType
  import args 1 bindings.GraphSignal.FftSection.windowSize
  import args 1 bindings.GraphSignal.FftSection.windowPaddSize
  import args 1 bindings.GraphSignal.FftSection.dynamic

  # Get the right time index
  x0    = sig.x0
  dx    = sig.dx
  index = int((@x-x0)/dx+.5)

  # Should test whether ..select1 exists, and return with no action if not
  l = [msge ${@objname}.^ list 'select*']
  if (l.length == 0) {
    printf "You should select a region first to specify the window size\n"
    return
  }
  s = l[0]
  {i j di dj} = [setg ..$s -rect]

  windowType     ='hamming'
  windowSize     = int(di/dx)
  windowSize = 2^(ceil(log(windowSize)/log(2)))
  windowPaddSize = min(2048,2*windowSize)
  dynamic        = 90

  # Get the spectrum
  window   =[new &signal]
  fftR     =[new &signal]
  fftI     =[new &signal]
  fft      =[new &signal]
  threshold=[new &signal]

  # Extract and window the piece of signal
  extract=sig[*b0,(index-windowSize//2):(index-windowSize//2+windowSize-1)]
  window = [stft window windowType windowSize]
  extract*=window

  # Zero padding
  extract1=extract[*b0,0:!windowPaddSize]

  # Fft
  fft extract1 fftR fftI
  fft=fftR*fftR+fftI*fftI

  # Size the dynamic range
  threshold = fft>= 10^(-dynamic/10)
  fft       = threshold*fft+10^(-dynamic/10)*(1-threshold)
  fft       = 10*log(fft)/log(10)

  # Then just display the energy
  xLabel = "Hz"
  yLabel = "dB"

  # Case the window already exists  
  if ([msge fftWindow exist]) {
    {x0 x1 y0 y1} = [setg fftWindow.fv1.view -bound]
    {imin imax}   = [stats minmax fft]
    ymin = fft.Y[imin]
    ymax = fft.Y[imax]
    update=0
    if (y0 > ymin) {
	y0=ymin
	update=1
    }
    if (y1 < ymax) {
	y1=ymax
	update=1
    }
    if (update) {
	setgu fftWindow.fv1 -graph1 fft  "1" -bound '*' '*' y0 y1
    } else {
	setgu fftWindow.fv1 -graph1 fft "1"
    }
  } else {
      disp fftWindow fft -..fv1 -yLabel yLabel -xLabel xLabel -..fv1.view -bound '*' '*' '?' '*'
# If we replace the above line with the following lines
# we are supposed to shift the fftWindow so as
# not to overlap the signal window.
# However there is a bug so far that locks the mouse
# in the signal window whenever the fftWindow is deleted

#    {hsize vsize}=[setg . -size]
#    {hpos vpos}  = [setg . -pos]
#      disp fftWindow fft -pos [= $hpos+$hsize] $vpos -..fv1 -yLabel yLabel -xLabel xLabel -..fv1.view -bound * * ? *
  }  
}

setbinding 'fftsection' "{{Shift + Right button = FFT spectrum}}"
setbinding 'fftsection' GraphSignal rightButtonDown shift {_FFTSection}
setbinding 'fftsection' GraphSignal rightButtonMotion shift {_FFTSection}

binding activate 'fftsection'