This notebook demonstarates the S-Method.

See Boashash+15 Section 6.2 for instance. last update: 9/30 (2017)



In [1]:

    
VERSION









    Out[1]:





v"0.6.1"



In [2]:

    
using PyPlot
import DSP



In [3]:

    
include("../juwvid.jl")









    Out[3]:





juwvid



In [4]:

    
# multicomponent data (Boashash+15,p346,Example 6.2.2)
nsample=1024
t,x=sampledata.genmultifm622(nsample);
fig=PyPlot.figure(figsize=(10,3))
PyPlot.plot(t,x)
PyPlot.xlabel("t")
PyPlot.ylabel("y")
PyPlot.savefig("triple.png")
#writedlm("mf622.txt",[t x])



In [5]:

    
# STFT
tfrstft=stft.tfrstft(x,NaN,NaN,NaN,NaN,NaN,4);









    



Use fft.



In [6]:

    
# pseudo Wigner Ville
z=DSP.Util.hilbert(x);
tfr=cohenclass.tfrpwv(z);









    



Single pseudo Wigner Ville
Use fft.



In [7]:

    
# S-method (constant Lp=6)
sm1=smethod.tfrsm(x,NaN,6,NaN,4);









    



Single S-method
Use fft.



In [8]:

    
#frequency width
2*6/nsample/(t[2]-t[1])/2









    Out[8]:




2.9970703125000027



In [9]:

    
# S-method (constant Lp=15)
sm2=smethod.tfrsm(x,NaN,15,NaN,4);









    



Single S-method
Use fft.

STFT is not affected by the cross terms, but, the resolution is poor. For the multi component FM, (pseudo) Wigner Ville has severe contaminations from cross terms. S-method has both advantages of the WV and STFT



In [10]:

    
fig=PyPlot.figure()
ax = fig[:add_subplot](2,2,1)
a=juwplot.tfrshow(abs.(tfr),t[2]-t[1],t[1],t[end],NaN,NaN,0.7*2,"CMRmap")
PyPlot.title("pseudo Wigner Ville")
PyPlot.ylabel("frequency")
PyPlot.ylim(0,120)
ax = fig[:add_subplot](2,2,2)
a=juwplot.wtfrshow(abs.(sm2),t[2]-t[1],t[1],t[end],NaN,NaN,0.7*2,"CMRmap")
PyPlot.title("S-method (Lp=15)")
PyPlot.ylim(0,120)
ax = fig[:add_subplot](2,2,3)
a=juwplot.wtfrshow(abs.(sm1),t[2]-t[1],t[1],t[end],NaN,NaN,0.7*2,"CMRmap")
PyPlot.title("S-method (Lp=6)")
PyPlot.xlabel("time")
PyPlot.ylabel("frequency")
PyPlot.ylim(0,120)
ax = fig[:add_subplot](2,2,4)
a=juwplot.wtfrshow(abs.(tfrstft).^2,t[2]-t[1],t[1],t[end],NaN,NaN,0.7*2,"CMRmap")
PyPlot.title("Spectrogram")
PyPlot.xlabel("time")
PyPlot.ylim(0,120)









    












    Out[10]:





(0, 120)

SM with N/2 = WV



In [11]:

    
# Alias-free PWV = SM for Lp=N/2
trfafwv=smethod.tfrsm(x,NaN,floor(Int,nsample/2),NaN,4);









    



Single S-method
Use fft.



In [12]:

    
a=juwplot.wtfrshow(abs.(trfafwv),t[2]-t[1],t[1],t[end],NaN,NaN,0.7*2,"CMRmap",0,600)
PyPlot.colorbar(a,shrink=0.7)
PyPlot.ylim(0,120)









    












    Out[12]:





(0, 120)



In [ ]: