;program to create nice plot as a GIF-image to be used for SO-presentation
pro mp,x,f,color=color,name=name,sz=sz,_extra=_extra, $
       linestyle=linestyle,thick=thick,noerase=noerase, $
       open=open,close=close,nodata=nodata,psym=psym
  common filename,file
  
  if keyword_set(close) eq 0 then begin
    
    if n_elements(xrange) eq 0 then xrange=[min(x),max(x)]
    if n_elements(yrange) eq 0 then yrange=[min(f),max(f)]
    if n_elements(sz) eq 0 then sz=[800,600]
    
    if keyword_set(noerase) eq 0  then begin
      file='../figures/'+name+'.png'
      window,xsize=sz(0),ysize=sz(1)
      
      userlct,/nologo
      !p.font=1
    endif
    !p.charsize=1
    
    plot,_extra=_extra, $
      xrange=xrange,yrange=yrange,/nodata,color=!p.color, $
      [xrange(0),xrange(1)],[yrange(0),yrange(1)]
    
    szf=size(f)
    if szf(0) eq 1 then np=1 else np=szf(2)
    if n_elements(color) ne np then color=9+intarr(np)
    if n_elements(linestyle) ne np then linestyle=intarr(np)
    if n_elements(psym) ne np then psym=intarr(np)
    if n_elements(thick) ne np then thick=intarr(np)+1
    for i=0,np-1 do begin
      if keyword_set(nodata) eq 0 then $
        oplot,x,f(*,i),color=color(i),_extra=_extra, $
        thick=thick(i),linestyle=linestyle(i),psym=psym(i)
    endfor
  endif
  
  if keyword_set(open) eq 0 or keyword_set(close) then begin
    img=tvrd(true=1)
;  tvlct,r,g,b,/get
    write_png,file,img
    print,'Image: '+file
  endif
  
end

function noise,x,rg
  
  noise=randomu(seed,n_elements(x))
  noise=noise*(max(rg)-min(rg))+min(rg)
  
  return,noise
end

pro sinusoidal
  
  nel=1000
  
  x=findgen(nel)/(nel-1)*2*!pi*5
  y1=sin(x)
  y2=x*0+1
  y3=x*0-1
  f=[[y1],[y2],[y3]]
  color=[1,9,9]
  
  mp,x,f,color=color,name='sinusoidal',yrange=[-1.1,1.1],/xst,/yst, $
    linestyle=[0,2,2],thick=[2,1,1],charsize=2,position=[.20,.25,.95,.95], $
    xtitle='Time',ytitle='Amplitude',sz=[550,350]
  
  x=findgen(nel)/(nel+1)*2
  y1=x*0 & y1(nel/2)=1.
  y2=x*0+1
  f=[[y1],[y2]]
  color=[1,9]
  mp,x,f,color=color,name='sinusoidal_f',yrange=[0,1.1],/xst,/yst, $
    linestyle=[0,2],thick=[2,1],charsize=2,position=[.22,.25,.95,.95], $
    xtitle='Frequency',ytitle='Amplitude',sz=[400,350],xtickname=[' ',' ','f0',' ']
  
  
end

pro complex_periodic
  
  nel=1000
  
  x=findgen(nel)/(nel-1)*2*!pi/2.
  
  a=[1.,.1,1.5,0.6]
  b=[0.3,1,0.5,1.6]
  n=findgen(n_elements(a))+1
  f1=1
  a0=0.5
  y1=a0/2. + (a##cos(2*!pi*f1*x#n)+b##sin(2*!pi*f1*x#n))
  
  y2=x*0
  f=[[y1],[y2]]
  color=[1,9]
  
  mp,x,f,color=color,name='complex_periodic',/xst, $
    linestyle=[0,1],thick=[2,1],charsize=2,position=[.20,.25,.95,.95], $
    xtitle='Time',ytitle='Amplitude',sz=[550,350]
  wait,1
  
  x=findgen(nel)/(nel-1)*(n_elements(a)+1)
  y=x*0
  ab=sqrt(a^2+b^2)
  for i=0,n_elements(a)-1 do begin
    dummy=min(abs(x-n(i))) & y(!c)=ab(i)
  endfor
  mp,x,y,color=1,name='complex_periodic_f',/xst, $
    linestyle=[0,1],thick=[2,1],charsize=2,position=[.22,.25,.95,.95], $
    xtitle='Frequency',ytitle='Amplitude',sz=[400,350], $
    xtickname=['0','1','2','3','4','5']+'f!L1!N'

end

pro almost_periodic
  
  nel=1000
  
  x=findgen(nel)/(nel-1)*2*!pi/2
  x0=[3,2,4]
  f0=[2,3,sqrt(50.)]
  th0=[30.,67.,128.]*!dtor
  y1=x*0
  for i=0,n_elements(x0)-1 do $
    y1=y1+ x0(i) * sin(2*!pi*f0(i)*x+th0(i))
  
  y2=x*0
  f=[[y1],[y2]]
  color=[1,9]
  
  mp,x,f,color=color,name='almost_periodic',/xst, $
    linestyle=[0,1],thick=[2,1],charsize=2,position=[.20,.25,.95,.95], $
    xtitle='Time',ytitle='Amplitude',sz=[550,350]
  wait,1
  
  x=findgen(nel)/(nel-1)*(max(f0)+1)
  y=x*0
  for i=0,n_elements(f0)-1 do begin
    dummy=min(abs(x-f0(i))) & y(!c)=x0(i)
  endfor

  mp,x,y,color=1,name='almost_periodic_f',/xst, $
    linestyle=[0,1],thick=[2,1],charsize=2,position=[.22,.25,.95,.95], $
    xtitle='Frequency',ytitle='Amplitude',sz=[400,350], $
    xtickname='f!L'+n2s(indgen(n_elements(f0))+1)+'!N', $
    xtickv=f0,xticks=n_elements(f0)
end

pro transient
  
  chsz=1.5
  
  nel=100000
  T=50
  dt=2*!pi*T/nel
  
  x=findgen(nel)/(nel-1)*2*!pi*T
  f=findgen(nel)/nel/dt
  
  
  n21=nel/2+1
  f=findgen(nel)
  f(n21)=n21-nel+findgen(n21-2)
  f=f/(nel*dt)
  
  y=5*exp(-.1*x)
;  y=cos(x*2*!pi*.5)
  mp,x,y,color=1,name='transient_exp',/xst,yrange=[0,5.1],$
    linestyle=[0,1],thick=[2,1],charsize=chsz,position=[.22,.25,.95,.95], $
    xtitle='Time',ytitle='Amplitude',sz=[300,200],xrange=[0,T]
  wait,1
  fy=sqrt((fft(y,-1))^2)
  mp,shift(f,-n21),shift(fy,-n21),color=1,name='transient_exp_f',/xst, $
    linestyle=[0,1],thick=[2,1],charsize=chsz,position=[.22,.25,.95,.95], $
    xtitle='Frequency',ytitle='Amplitude',sz=[300,200],xrange=[0,0.1]
  wait,1
  
  y=5*exp(-.1*x)*cos(2*!pi*.1*x)
  mp,x,y,color=1,name='transient_expcos',/xst,yrange=[-3,5.1],$
    linestyle=[0,1],thick=[2,1],charsize=chsz,position=[.22,.25,.95,.95], $
    xtitle='Time',ytitle='Amplitude',sz=[300,200],xrange=[0,T]
  wait,1
  fy=sqrt((fft(y,-1))^2)
  mp,shift(f,-n21),shift(fy,-n21),color=1,name='transient_expcos_f',/xst, $
    linestyle=[0,1],thick=[2,1],charsize=chsz,position=[.22,.25,.95,.95], $
    xtitle='Frequency',ytitle='Amplitude',sz=[300,200],xrange=[0,0.5]
  wait,1
  
  y=x*0
  y(0:nel/30)=1.
  mp,x,y,color=1,name='transient_box',/xst,yrange=[0,1.1],$
    linestyle=[0,1],thick=[2,1],charsize=chsz,position=[.22,.25,.95,.95], $
    xtitle='Time',ytitle='Amplitude',sz=[300,200],xrange=[0,T]
  wait,1
  fy=sqrt((fft(y,-1))^2)
  mp,shift(f,-n21),shift(fy,-n21),color=1,name='transient_box_f',/xst, $
    linestyle=[0,1],thick=[2,1],charsize=chsz,position=[.22,.25,.95,.95], $
    xtitle='Frequency',ytitle='Amplitude',sz=[300,200],xrange=[0,0.3]
  wait,1
end

pro random
  common greek_letters
  
  greek
  nel=200
  T=10.
  x=findgen(nel)/(nel-1)*T
  chsz=0.5
  sm=9
  
  y0=x*0.
  
  !p.multi=[0,1,3]
  
  for i=0,2 do begin
    ue=(lindgen(nel))(sm:nel-sm)
    y1=smooth(randomu(seed,nel)-0.3,9)
    mp,x(ue),[[y0(ue)],[y1(ue)]], $
      color=[9,1],name='random',/xst,yrange=[-0.1,.6],$
      linestyle=[1,0],thick=[1,2],charsize=chsz, $
      xtitle=' ',ytitle=' ',sz=[400,600],xrange=[0,T], $
      xtickname=strarr(32)+' ',ytickname=strarr(32)+' ',noerase=i gt 0,/open
    t1=nel/3. & t2=nel/3.*2
    oplot,color=9,linestyle=2,x(t1)+[0,0],[!y.crange(0),y1(t1)]
    oplot,color=9,linestyle=2,x(t2)+[0,0],[!y.crange(0),y1(t2)]
    xyouts,/data,alignment=0.5,x(t1),!y.crange(1),'!C!Ct!L1!N',charsize=2.
    xyouts,/data,alignment=0.5,x(t2),!y.crange(1),'!C!Ct!L1!N + '+f_tau,charsize=2.
  endfor
  mp,/close
  
  !p.multi=0
end

pro cc_timedelay
  
  !p.multi=0
  chsz=2
  nel=200
  T=10.
  x=findgen(nel)/(nel-1)*T
  xi=indgen(nel)
  t0=0.5
  x1=3.
  x2=7.
  
  
  
  signal1=(1.*exp(-(x-x1)^2/t0) + noise(x,[-.1,+.1]) +1)
  signal2=(1.*exp(-(x-x2)^2/t0) + noise(x,[-.5,+.1]) +0)
  
  mp,x,[[signal1],[signal2]],color=[1,2], $
    name='cc_timedelay_signal',/xst,$
    linestyle=[0,0],thick=[2,2],charsize=chsz,position=[.15,.15,.95,.95], $
    xtitle='Time',ytitle='Signal',sz=[550,400],xrange=[0,T]
  
  lag=lindgen(nel)*2-2
  cc=fltarr(nel*2-3)
  cc1=fltarr(nel*2-3)
  xcc=fltarr(nel*2-3)
  for lg=-(nel-2),nel-2 do begin
    cc(lg+nel-2)=c_correlate(signal1,signal2,lg)
    xcc(lg+nel-2)=x(abs(lg))*([1,-1])(lg lt 0)
    
    olap=where(xi+lg ge 0 and xi+lg lt nel)
    if lg le 0 then begin
      s1=signal1(indgen(n_elements(olap)))
      s2=signal2(olap)
    endif else begin
      s1=signal1(lg+indgen(n_elements(olap)))
      s2=signal2(indgen(n_elements(olap)))
    endelse
    cc1(lg+nel-2)=1./n_elements(olap)*total(s1*s2) 
    
;    plot,xcc(indgen(n_elements(olap))),s1 & oplot,color=1,xcc(indgen(n_elements(olap))),s2; & wait,.1
;print,lg+nel-2,lg,n_elements(olap),cc1(lg+nel-2)

  endfor
; cc=cc1
  
  wait,1
  mp,xcc,cc,color=1, $
    name='cc_timedelay_sig1sig2',/xst,$
    linestyle=0,thick=2,charsize=chsz,position=[.15,.15,.95,.95], $
    xtitle='Time',ytitle='Cross Correlation',sz=[550,400],/open
  
  dummy=max(cc) & im=!c
  
  oplot,color=9,linestyle=2,xcc(im)+[0,0],[!y.crange(0),cc(im)]
  
  xyouts,/data,xcc(im),!y.crange(1),charsize=2.,alignment=0.5, $
    '!Clag = '+string(xcc(im),format='(f5.2)')
  
  mp,/close
  
                                ;use fft for computation
  fft1=fft(signal1)
  fft2=fft(signal2)
  
  fft_cc=fft2*conj(fft1)
  cc=(fft(fft_cc,1))
  
  wait,1
  mp,x,cc,color=1, $
    name='cc_fft_timedelay_sig1sig2',/xst,$
    linestyle=0,thick=2,charsize=chsz,position=[.15,.15,.95,.95], $
    xtitle='Time',ytitle='Cross Correlation',sz=[550,400],/open
  
  dummy=max(cc) & im=!c
  
  oplot,color=9,linestyle=2,x(im)+[0,0],[!y.crange(0),cc(im)]
  
  xyouts,/data,x(im),!y.crange(1),charsize=2.,alignment=0.5, $
    '!Clag = '+string(x(im),format='(f5.2)')
  
  mp,/close
  
end

pro cc_signoise
  
  !p.multi=0
  chsz=2
  nel=200
  T=10.
  x=findgen(nel)/(nel-1)*T
  xi=indgen(nel)
  t0=0.5
  x1=T/3.
  x2=T/2.
  x3=T/1.5
  A0=0.1

  signal1=(5*A0*exp(-(x-x1)^2/t0) + noise(x,[-.0,+.0]) +0)+ $
    (3*A0*exp(-(x-x2)^2/t0) + noise(x,[-.0,+.0]) +0) + $
    (1*A0*exp(-(x-x3)^2/t0) + noise(x,[-.0,+.0]) +0)
; signal2=(A0*exp(-(x-x2)^2/t0) + noise(x,[-.2,+.2]) +0)
  signal2=signal1/2. + noise(x,[-.1,+.1]) 
  signal3=signal1/5. + noise(x,[-.1,+.1]) 
  
  mp,x,[[signal1],[signal2],[signal3]],color=[1,2,3], $
    name='cc_signoise_signal',/xst,$
    linestyle=[0,0,0],thick=[2,2,2],charsize=chsz,position=[.15,.15,.95,.95], $
    xtitle='Time',ytitle='Signal',sz=[550,400],xrange=[0,T]
  
  
  lag=lindgen(nel)*2-2
  cc1=fltarr(nel*2-3)
  cc2=fltarr(nel*2-3)
  xcc=fltarr(nel*2-3)
  for lg=-(nel-2),nel-2 do begin
    cc1(lg+nel-2)=c_correlate(signal1,signal2,lg)
    cc2(lg+nel-2)=c_correlate(signal1,signal3,lg)
    xcc(lg+nel-2)=x(abs(lg))*([1,-1])(lg lt 0)
  endfor
  
  wait,1
  mp,xcc,[[cc1],[cc2]],color=[2,3], $
    name='cc_signoise_sig1sig2',/xst,$
    linestyle=[0,0],thick=[2,2],charsize=chsz,position=[.15,.15,.95,.95], $
    xtitle='Time',ytitle='Cross Correlation',sz=[550,400],/open
  
  dummy=max(cc1) & im=!c
  oplot,color=9,linestyle=2,xcc(im)+[0,0],[!y.crange(0),cc1(im)]
  xyouts,/data,xcc(im),!y.crange(1),charsize=2.,alignment=0.5, $
    '!Clag = '+string(xcc(im),format='(f5.2)')
  
  mp,/close
end


pro dft
  
  chsz=2
  n2=500
  N=2*n2
  order=n2                ;use FFT reconstruction only until order=...
  T=2*!pi*2.5
  x=findgen(n)/(n-1)*T
  y=x*0
  y(where(x ge 2 and x le 6))=1.5
  y(where(x ge 3 and x le 5))=2.
y=1*sin(5*x)
  y=y
  
                                ;calculate dft coefficients
  t0=systime(1)
  aq=fltarr(n2+1)
  bq=fltarr(n2)
  aq(0)=1./n*total(y)
  finn=(findgen(n)+1)
  aq(n2)=1./n*total(y*cos(finn*!pi))
  for q=1,n2-1 do begin
    al=2*!pi*q*finn/N
    aq(q)=2./N*total(y*cos(al))
    bq(q)=2./N*total(y*sin(al))
  endfor
  
                                ;reconstruct y from fourier coefficients
  df=x*0.
  order=order<n2
  finq=(findgen(order+1))
  for k=1,N do begin
    al=2*!pi*finq*k/N
    df(k-1)=aq(0)+(total(aq(1:order  )*cos(al(1:order  )))+ $
                   total(bq(1:order-1)*sin(al(1:order-1))) )
  endfor
  t_dft=systime(1)-t0
  
                                ;run fft
  t0=systime(1)
  ffy=fft(y)
  ff_full=fft(ffy,1)
  t_fft=systime(1)-t0
  
                                ;reconstruct fft only to specified
                                ;order
  ff=x*0.
  finq=findgen(order)
  for k=0,n-1 do begin
    ff(k)=(total(ffy(finq)*exp(complex(0,1)*2*!pi*finq*k/n)) +  $
           total(ffy(N-1-finq)*exp(complex(0,1)*2*!pi*(N-1-finq)*k/n)))
  endfor
  
  print,t_dft,t_fft,t_dft/t_fft
  xm=!x.margin
  ym=!y.margin
  !p.multi=[0,1,2]
  !x.margin=xm/(chsz*.9)
  !y.margin=ym/(chsz*.7)
  mp,x,[[y],[df],[ff]],color=[1,2,3], $
    name='dft',/xst,$
    linestyle=[0,0,0],thick=[2,2,2],charsize=chsz, $
    xtitle='Time',ytitle='Signal',sz=[600,500],/open
  
  dx=x(1)-x(0)
  n21=n/2+1
  f=findgen(n)
  f[n21]=n21-n+findgen(n21-2)
  f=f/n/dx
  
  sf=shift(f,-n21)
  sffy=shift(ffy,-n21)
  
  mp,sf,abs([[sffy],[sffy]]),color=[1,2], $
    name='dft',/xst, $
    xrange=[0,f(order)]*.1,$
    linestyle=[0,0],thick=[2,2],charsize=chsz, $
    xtitle='Frequency',ytitle='Signal',sz=[600,500],/noerase,/open,psym=4
  for i=0,n_elements(sf)-1 do $
    oplot,color=1,sf(i)+[0,0],[0,abs(sffy(i))],linestyle=0
  
  dummy=max(abs(sffy(n21:*))) & pm=!c+n21
  oplot,color=9,linestyle=2,sf(pm)+[0,0],[0,abs(sffy(pm))]
  xyouts,/data,sf(pm),!y.crange(1),alignment=0.5,charsize=chsz, $
    '!C!Cf='+string(sf(pm),format='(f5.2)')+ $
    '!Ct='+string(1./sf(pm),format='(f6.2)')
  
  mp,/close
  
  
  !p.multi=0
  !x.margin=xm
  !y.margin=ym
  stop
end

pro spec_leakage
  
  !x.margin=[10,3] & !y.margin=[4,2] 

  chsz=2
  n2=500
  N=2*n2
  
  for il=0,2 do begin
    if il eq 0 then T=4 $
    else            T=4.5
  x=findgen(n)/(n-1)*T
  
  ysig=1*sin(x*2*!pi)
  w=exp(-0.5*(.8*(x-T/2.))^2)
  if il eq 2 then y=ysig*w $
  else y=ysig
  
  
                                ;run fft
  ffy=fft(y)
  ff=fft(ffy,1)
  
  xm=!x.margin
  ym=!y.margin
  !p.multi=[0,1,2]
  !x.margin=xm/(chsz*.9)
  !y.margin=ym/(chsz*.7)
  mp,x,y,color=[1], $
    name='spec_leakage'+n2s(il),/xst,$
    yrange=[-1.1,1.1], $
    linestyle=[0],thick=[2],charsize=chsz, $
    xtitle='Time',ytitle='Signal',sz=[500,500],/open
  
  if il eq 2 then begin
    oplot,color=2,x,ysig
    oplot,color=3,x,w
  endif
  
  dx=x(1)-x(0)
  n21=n/2+1
  f=findgen(n)
  f[n21]=n21-n+findgen(n21-2)
  f=f/n/dx
  
  sf=shift(f,-n21)
  sffy=shift(ffy,-n21)
  
  mp,sf,[[sf*0],[abs(sffy)]],color=[9,1],/xst, $
    xrange=[0,max(f)]*.02,$
    yrange=max(abs(sffy))*.1*[-1,1]+[0,max(abs(sffy))], $
    linestyle=[1,0],thick=[1,2],charsize=chsz, $
    xtitle='Frequency',ytitle='Signal',/noerase,/open,psym=[-0,4]
  for i=0,n_elements(sf)-1 do $
    oplot,color=1,sf(i)+[0,0],[0,abs(sffy(i))],linestyle=0
  
  
  mp,/close
  
  wait,.5
  
  !p.multi=0
  !x.margin=xm
  !y.margin=ym
  endfor
end

pro ipol
  
  chsz=2
  nel =100
  T=4.
  x=findgen(nel)/(nel-1)*T
  y=sin(2*!pi*x)+noise(x,[-1,1]*.2)
  
  sz=[500,400]
  ps=[.15,.15,.95,.95]
  
  mp,x,[[y*0],[y]],color=[9,1],name='interpol1',/xst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz,position=ps, $
    xtitle='Time',ytitle='Amplitude',sz=sz,yrange=[-1.25,1.25],/yst, $
    xrange=[0,4]
  
  ffy=fft(y)
  dx=x(1)-x(0)
  n21=nel/2+1
  f=findgen(nel)
  f[n21]=n21-nel+findgen(n21-2)
  f=f/nel/dx
  sf=shift(f,-n21)
  sffy=shift(ffy,-n21)
  
  wait,.5
  mp,sf,[[sf*0],[abs(sffy)]],color=[9,1],/xst,position=[.15,.20,.95,.95], $
    name='interpol1_ff',sz=[sz(0),.8*sz(1)], $
    xrange=[0,1/T*20],$
;    yrange=max(abs(sffy))*.1*[-1,1]+[0,max(abs(sffy))], $
    yrange=[-.1,.7], $
    linestyle=[1,0],thick=[1,2],charsize=chsz, $
    xtitle='Frequency',ytitle='Signal',/open,psym=[-0,4]
  for i=0,n_elements(sf)-1 do $
    oplot,color=1,sf(i)+[0,0],[0,abs(sffy(i))],linestyle=0
  mp,/close
  
  ns=25
  spx=(findgen(ns)+randomu(seed,ns))/ns*(max(x)-min(x)) + min(x)
  idx=lonarr(ns)
  for i=0,ns-1 do begin
    dummy=min(abs(x-spx(i)))
    idx(i)=!c
  endfor
    
  wait,.5
  mp,x,[[y*0],[y]],color=[9,1],name='interpol2',/xst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz,position=ps, $
    xtitle='Time',ytitle='Amplitude',sz=sz,/open,yrange=[-1.25,1.25],/yst, $
    xrange=[0,4]
  
  oplot,color=3,psym=4,spx,y(idx)
  for i=0,ns-1 do $
    oplot,color=3,linestyle=0,[0,0]+spx(i),[0,y(idx(i))]
  
  mp,/close
  
  wait,.5
  mp,x,[[y*0],[y]],color=[9,1],name='interpol3',/xst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz,position=ps, $
    xtitle='Time',ytitle='Amplitude',sz=sz,/open,/nodata, $
    yrange=[-1.25,1.25],/yst, $
    xrange=[0,4]
  
  oplot,color=9,linestyle=1,!x.crange,[0,0]
  oplot,color=3,psym=4,spx,y(idx)
  for i=0,ns-1 do $
    oplot,color=3,linestyle=0,[0,0]+spx(i),[0,y(idx(i))]
  
  mp,/close
  
  dx=abs(spx-shift(spx,1))
;  mindx=min(dx)
  mindx=total((dx(sort(dx)))(1:5))/5.
  nn=(max(spx)-min(spx))/mindx
  xn=findgen(nn)/(nn-1)*(max(spx)-min(spx)) + min(spx)
  
  add=['linear','lsquadratic','spline','quadratic']
  for ii=0,3 do begin
    yi=interpol(y(idx),spx,xn,quadratic=ii eq 3,lsq=ii eq 1,spline=ii eq 2)
  
  wait,.5
  mp,xn,[[yi*0],[yi]],color=[9,2],name='interpol4'+add(ii),/xst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz,position=ps, $
    xtitle='Time',ytitle='Amplitude',sz=sz,/open,yrange=[-1.25,1.25],/yst, $
    xrange=[0,4]
  oplot,color=3,psym=4,spx,y(idx)
  oplot,color=1,linestyle=0,x,y
  for i=0,ns-1 do $
    oplot,color=3,linestyle=0,[0,0]+spx(i),[0,y(idx(i))]
  mp,/close
endfor
  
  wait,.5
  mp,xn,[[yi*0],[yi]],color=[9,2],name='interpol5',/xst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz,position=ps, $
    xtitle='Time',ytitle='Amplitude',sz=sz,/open,yrange=[-1.25,1.25],/yst, $
    xrange=[0,4]
  oplot,color=3,psym=4,xn,yi
  for i=0,nn-1 do $
    oplot,color=3,linestyle=0,[0,0]+xn(i),[0,yi(i)]
  mp,/close
  
  ffy=fft(yi)
  dx=xn(1)-xn(0)
  n21=nn/2+1
  f=findgen(nn)
  f[n21]=n21-nn+findgen(n21-2)
  f=f/nn/dx
  sf=shift(f,-n21)
  sffy=shift(ffy,-n21)
  
  wait,.5
  mp,sf,[[sf*0],[abs(sffy)]],color=[9,1],/xst,position=[.15,.20,.95,.95], $
    name='interpol6_ff',sz=[sz(0),.8*sz(1)], $
    xrange=[0,1/T*20],$
;    yrange=max(abs(sffy))*.1*[-1,1]+[0,max(abs(sffy))], $
    yrange=[-.1,.7], $
    linestyle=[1,0],thick=[1,2],charsize=chsz, $
    xtitle='Frequency',ytitle='Signal',/open,psym=[-0,4]
  for i=0,n_elements(sf)-1 do $
    oplot,color=1,sf(i)+[0,0],[0,abs(sffy(i))],linestyle=0
  mp,/close
   
end

pro denoise
  
  chsz=2
  nel =150
  T=8.
  x=findgen(nel)/(nel-1)*T
  y=.5*sin(2*!pi*x)+.7*sin(3*!pi*x)+noise(x,[-1,1]*.5)
  
  sz=[500,400]
  ps=[.15,.15,.95,.95]
  
  mp,x,[[y*0],[y]],color=[9,1],name='denoise1',/xst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz,position=ps, $
    xtitle='Time',ytitle='Amplitude',sz=sz,yrange=[-1.5,1.5],/yst, $
    xrange=[0,T]
  
  ffy=fft(y)
  dx=x(1)-x(0)
  n21=nel/2+1
  f=findgen(nel)
  f[n21]=n21-nel+findgen(n21-2)
  f=f/nel/dx
  sf=shift(f,-n21)
  sffy=shift(ffy,-n21)
  
  wait,.5
  mp,sf,[[sf*0],[abs(sffy)]],color=[9,1],/xst,position=ps, $
    name='denoise1_ff',sz=sz, $
    xrange=[0,1/T*nel/2.],$
    yrange=[-.1,.5],/yst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz, $
    xtitle='Frequency',ytitle='Signal',/open,psym=[-0,4]
 for i=0,n_elements(sf)-1 do $
    oplot,color=1,sf(i)+[0,0],[0,abs(sffy(i))],linestyle=0
  mp,/close
  
  cutoff=1./t*nel/2. * 0.4
  cidx=where(f ge cutoff or f le -cutoff)
  ffy(cidx)=0.
  sffy=shift(ffy,-n21)

  wait,.5
  mp,sf,[[sf*0],[abs(sffy)]],color=[9,1],/xst,position=ps, $
    name='denoise2_ff',sz=sz, $
    xrange=[0,1/T*nel/2.],$
    yrange=[-.1,.5],/yst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz, $
    xtitle='Frequency',ytitle='Signal',/open,psym=[-0,4]
  for i=0,n_elements(sf)-1 do $
    oplot,color=1,sf(i)+[0,0],[0,abs(sffy(i))],linestyle=0
  oplot,color=9,linestyle=2,cutoff+[0,0],!y.crange
  mp,/close
 
  ff=fft(ffy,1)
  
  wait,.5
  mp,x,[[y*0],[y],[ff]],color=[9,8,1],name='denoise2',/xst, $
    linestyle=[1,0,0],thick=[1,1,2],charsize=chsz,position=ps, $
    xtitle='Time',ytitle='Amplitude',sz=sz,yrange=[-1.5,1.5],/yst, $
    xrange=[0,T]
 
  cutoff=1./t*nel/2. * 0.4
  threshold=0.1
  cidx=where(ffy le threshold)
  ffy(cidx)=0.
  sffy=shift(ffy,-n21)

  wait,.5
  mp,sf,[[sf*0],[abs(sffy)]],color=[9,1],/xst,position=ps, $
    name='denoise3_ff',sz=sz, $
    xrange=[0,1/T*nel/2.],$
    yrange=[-.1,.5],/yst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz, $
    xtitle='Frequency',ytitle='Signal',/open,psym=[-0,4]
  for i=0,n_elements(sf)-1 do $
    oplot,color=1,sf(i)+[0,0],[0,abs(sffy(i))],linestyle=0
  oplot,color=9,linestyle=2,!x.crange,[0,0]+threshold
  mp,/close
 
  ff=fft(ffy,1)
  
  wait,.5
  mp,x,[[y*0],[y],[ff]],color=[9,8,1],name='denoise3',/xst, $
    linestyle=[1,0,0],thick=[1,1,2],charsize=chsz,position=ps, $
    xtitle='Time',ytitle='Amplitude',sz=sz,yrange=[-1.5,1.5],/yst, $
    xrange=[0,T]
 
end

pro chirp
  
  sz=[500,400]
  ps=[.15,.15,.95,.95]
  
  chsz=2
  nel =500
  T=8.
  x=findgen(nel)/(nel-1)*T
  f=2*!pi*1.*x/3
  y=1*sin(f*x)
  
  
  mp,x,[[y*0],[y]],color=[9,1],name='chirp',/xst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz,position=ps, $
    xtitle='Time',ytitle='Amplitude',sz=sz,yrange=[-1.5,1.5],/yst, $
    xrange=[0,T]
  
  ffy=fft(y)
  dx=x(1)-x(0)
  n21=nel/2+1
  f=findgen(nel)
  f[n21]=n21-nel+findgen(n21-2)
  f=f/nel/dx
  sf=shift(f,-n21)
  sffy=shift(ffy,-n21)
  
  wait,.5
  mp,sf,[[sf*0],[abs(sffy)]],color=[9,1],/xst,position=ps, $
    name='chirp_ff',sz=sz, $
    xrange=[0,1/T*nel/5.],$
    yrange=[-.02,.12],/yst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz, $
    xtitle='Frequency',ytitle='Signal',/open,psym=[-0,4]
  for i=0,n_elements(sf)-1 do $
    oplot,color=1,sf(i)+[0,0],[0,abs(sffy(i))],linestyle=0
  mp,/close
  
  
end


pro wl_demo
  
  sz=[500,300]
  ps=[.15,.22,.95,.95]
  
  chsz=2
  nel =10000
  T=12.234345636
  x=findgen(nel)/(nel-1)*T
;  f=2*!pi*1.*x/3
  f0=[6.,4.,2.,1]*2*!pi
  dt=[0.,2,4,8]
;  f0=[100.,50.,25.,10.]*2*!pi
;  dt=[0,.3,.6,.8]*T
  dt=[dt,T+1]
  y=x*0
  for i=0,n_elements(f0)-1 do begin
    idx=where(x ge dt(i) and x lt  dt(i+1))
    y(idx)=1.*cos(f0(i)*x(idx))
  endfor
  y1=total(1.*cos(f0 # x),1)/n_elements(f0)
  
  mp,x,[[y*0],[y]],color=[9,1],name='wl_demo1',/xst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz,position=ps, $
    xtitle='Time',ytitle='Amplitude',sz=sz,yrange=[-1.5,1.5],/yst, $
    xrange=[0,T]

  wait,.5
  mp,x,[[y1*0],[y1]],color=[9,1],name='wl_demo2',/xst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz,position=ps, $
    xtitle='Time',ytitle='Amplitude',sz=sz,yrange=[-1.5,1.5],/yst, $
    xrange=[0,T]
  
;  y=y*(1-abs(x-(T/2.))/(T/2.))  ;windowing with triangle
  ffy=fft(y)
  ffy1=fft(y1)
  dx=x(1)-x(0)
  n21=nel/2+1
  f=findgen(nel)
  f[n21]=n21-nel+findgen(n21-2)
  f=f/nel/dx
  sf=shift(f,-n21)
  sffy=shift(ffy,-n21)
  sffy1=shift(ffy1,-n21)
  
  wait,.5
  mp,sf,[[sf*0],[abs(sffy)]],color=[9,1],/xst,position=ps, $
    name='wl_demo1_ff',sz=sz, $
    xrange=[0,max(f0)/2/!pi*1.5], $
;    yrange=[-.02,.12],/yst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz, $
    xtitle='Frequency',ytitle='Signal',/open,psym=[-0,4]
  for i=0l,n_elements(sf)-1 do $
    oplot,color=1,sf(i)+[0,0],[0,abs(sffy(i))],linestyle=0
  mp,/close
  
  wait,.5
  mp,sf,[[sf*0],[abs(sffy1)]],color=[9,1],/xst,position=ps, $
    name='wl_demo2_ff',sz=sz, $
    xrange=[0,max(f0)/2/!pi*1.5], $
;    yrange=[-.02,.12],/yst, $
    linestyle=[1,0],thick=[1,2],charsize=chsz, $
    xtitle='Frequency',ytitle='Signal',/open,psym=[-0,4]
  for i=0l,n_elements(sf)-1 do $
    oplot,color=1,sf(i)+[0,0],[0,abs(sffy1(i))],linestyle=0
  mp,/close
  
  for i=0,n_elements(f0)-1 do begin
    dw=2.
    w=x ge dt(i) and x lt dt(i)+dw
    wait,.5
    mp,x,[[w*0],[w]],color=[9,1],name='stft_w'+n2s(i),/xst, $
      linestyle=[1,0],thick=[1,2],charsize=chsz,position=ps, $
      xtitle='Time',ytitle='Amplitude',sz=sz,yrange=[-1.5,1.5],/yst, $
      xrange=[0,T]
    
    wait,.5
    mp,x,[[w*0],[y*w]],color=[9,1],name='stft_yw'+n2s(i),/xst, $
      linestyle=[1,0],thick=[1,2],charsize=chsz,position=ps, $
      xtitle='Time',ytitle='Amplitude',sz=sz,yrange=[-1.5,1.5],/yst, $
      xrange=[0,T]
    
  endfor
  
end