up [pdf]
from rsf.proj import *
from math import *
from rsf.prog import RSFROOT
from rsf.recipes.beg import server as private


#----------------------## parameters of synthetic signal

nt=501
dt = 0.008
nf = 257
ot=0

wf = 2*pi

def Grey(data,other):
    Result(data,
        '''
     grey  color=j screenratio=0.8 scalebar=n 
    barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
    label2=Time unit2=s label1=Frequency unit1=Hz n2tic=20 grid=n pclip=99.085 wanttitle=n %s
    ''' %(other))

nlevel0=0
nf0=247
df0=0.244141
###########################################################
# EXAMPLE 0000---cos+cos+cos+spike
###########################################################

Flow('s1-0',None,
     '''
     math n1=%d o1=0 d1=%g n2=1 output="cos(%g*(10*x1))+cos(%g*(20*x1))+cos(%g*(30*x1))" | 
     put label1=Time unit1=s label2=Amplitude 
     ''' % (nt,dt,wf,wf,wf))
Flow('add-0',None,'spike n1=%d o1=0 d1=%g n2=1 nsp=2 mag=10,10 k1=%d,%d' % (nt,dt,nt/2,nt/2+40))
Flow('s-0','s1-0 add-0','add ${SOURCES[1]} | noise seed=2012 var=%g' % (nlevel0))

Result('s-0','graph pad=n screenratio=0.8 crowd1=0.75 label2=Amplitude label1=Time unit1=s crowd2=0.3 wanttitle=n')

## TF 15
Flow('proj-0','s-0','timefreq rect=15 | window max2=60')
Result('proj-0',
       '''
       transp | scale axis=2 | grey wanttitle=n color=j screenratio=0.8 scalebar=n bartype=v
       barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
       label2=Time unit2=s label1=Frequency unit1=Hz n2tic=20 grid=n
       ''')
## TF 30
Flow('proj-0-30','s-0','timefreq rect=30 | window max2=60')
Result('proj-0-30',
       '''
       transp | scale axis=2 | grey wanttitle=n color=j screenratio=0.8 scalebar=n bartype=v
       barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
       label2=Time unit2=s label1=Frequency unit1=Hz n2tic=20 grid=n
       ''')


## STFT transform
Flow('st-0','s-0','st | math output="abs(input)" | real |  window max2=60')
Result('st-0',
       '''
       transp | scale axis=2 | grey wanttitle=n color=j screenratio=0.8 scalebar=n bartype=v
       barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
       label2=Time unit2=s label1=Frequency unit1=Hz n2tic=20 grid=n
       ''')

########################################################################
# INITIALIZATION for Matlab
########################################################################
matlab         = WhereIs('matlab')
matROOT = '../Matfun'
matfun = 'Synth'
matlabpath = os.environ.get('MATLABPATH',os.path.join(RSFROOT,'lib'))

if not matlab:
    sys.stderr.write('\nCannot find Matlab.\n')
    sys.exit(1)

fs=1/dt         #sampling frequency
nf=257          #number of frequency samples
df=fs/2/(nf-1)  #frequency interval

############################################################
## generate and process synthetic data
############################################################
Flow('tfsswt0 tfwt0',[os.path.join(matROOT,matfun+'.m'),'s-0'],
     '''MATLABPATH=%(matlabpath)s %(matlab)s 
     -nosplash -nojvm -r "addpath %(matROOT)s;%(matfun)s('${SOURCES[1]}','${TARGET}','${TARGETS[1]}',%(dt)g,%(nt)d,%(nf)d);quit"
     '''%vars(),stdin=0,stdout=-1)

Flow('gch1-tfsswt','tfsswt0','put o1=%g d1=%g n1=%d o2=%g d2=%g n2=%d'%(0,df,nf,0,dt,nt))
Flow('gch1-tfwt','tfwt0','put o1=%g d1=%g n1=%d o2=%g d2=%g n2=%d'%(0,df,nf,0,dt,nt))

Grey('gch1-tfsswt','')
Grey('gch1-tfwt','')

End()

sfmath
sfput
sfspike
sfadd
sfnoise
sfgraph
sftimefreq
sfwindow
sftransp
sfscale
sfgrey
sfst
sfreal