Comparison between the FFT mode and the non-unifrom FFT mode for the pseudo Wigner distribution



In [1]:

    
import DSP
using PyPlot



In [2]:

    
include("../juwvid.jl")









    Out[2]:





juwvid



In [3]:

    
## sin FM 
nsample=4096;
xs,ys,iws,ynorms=sampledata.genfm(nsample,2*pi,2*pi/365,1.0,365.0);
z=DSP.Util.hilbert(ys);
PyPlot.plot(xs,ys)









    












    Out[3]:





1-element Array{PyCall.PyObject,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7fa8e2eb0630>

the Pseudo Wigner Ville distribution with FFT



In [4]:

    
tfrpf=cohenclass.tfrpwv(z,NaN,NaN,NaN,NaN,NaN,0);









    



Single pseudo Wigner Ville
Use fft.



In [5]:

    
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.tfrshow(abs.(tfrpf),xs[2]-xs[1],xs[1],xs[end],NaN,NaN,0.7,"Spectral")
PyPlot.xlabel("time [day]")
PyPlot.ylabel("frequency [1/day]")
PyPlot.title("Pseudo WV")
ax = fig[:add_subplot](1,2,2)
a=juwplot.tfrshow(abs.(tfrpf),xs[2]-xs[1],xs[1],xs[end],NaN,NaN,0.7*130,"Spectral")
PyPlot.xlabel("time [day]")
PyPlot.title("Zoom-up")
PyPlot.ylim(0.98,1.02)









    












    Out[5]:





(0.98, 1.02)



In [6]:

    
# extracting the IF
indf=extif.maxif(abs.(tfrpf));
dx=xs[2]-xs[1]
fx=juwutils.index_to_frequency(indf,NaN,dx,nsample);









    



Assuming nft = nsample.



In [7]:

    
# We extract frequency indices of 0.98 and 1.02 [1/day] because we want to know the signal around 1. 
nnufft=1024
js,je=juwutils.frequency_to_index([0.98,1.02], xs[2]-xs[1], nsample,nnufft)









    Out[7]:





2-element Array{Float64,1}:
 715.575
 744.782

the Pseudo Wigner-Ville with the Zero-padding (usually slow)



In [8]:

    
Nz=100
tfrpfz=cohenclass.tfrpwv(z,NaN,NaN,[js,je],NaN,NaN,0,"zeropad",4,Nz);









    



Single pseudo Wigner Ville
Use Zero-padding fft.



In [9]:

    
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.tfrshow(abs.(tfrpf),xs[2]-xs[1],xs[1],xs[end],NaN,NaN,0.7*130,"Spectral")
PyPlot.xlabel("time [day]")
PyPlot.title("Zoom-up (FFT)")
PyPlot.ylim(0.985,1.02)
ax = fig[:add_subplot](1,2,2)
a=juwplot.tfrshow(abs.(tfrpfz),xs[2]-xs[1],xs[1],xs[end],round(Int,js),round(Int,je),0.7,"Spectral")
PyPlot.xlabel("time [day]")
PyPlot.title("zero-padding")
PyPlot.ylim(0.985,1.02)









    












    Out[9]:





(0.985, 1.02)



In [10]:

    
indfz=extif.maxif(abs.(tfrpfz));
fz=juwutils.index_to_frequency(indfz,NaN,xs[2]-xs[1],nsample,NaN,je,js,Nz);

PyPlot.title("no noise")
PyPlot.xlabel("time")
PyPlot.ylabel("instantaneous frequency")
PyPlot.ylim(0.995,1.005)
PyPlot.plot(xs,fx, color="gray",lw=3)
PyPlot.plot(xs,fz, color="C0",lw=3)
PyPlot.plot(xs,iws/(2*pi), color="C3", ls="dashed",lw=3)
PyPlot.legend(["IF from PWV (FFT)","IF from PWV (zeropad)","Input IF"])









    












    



Assuming nft = nsample.






    Out[10]:





PyObject <matplotlib.legend.Legend object at 0x7fa8daa445c0>

the Pseudo Wigner-Ville with NUFFT (faster than zero-padding)



In [11]:

    
fin=collect(linspace(js,je,nnufft));
tfrpfn=cohenclass.tfrpwv(z,NaN,NaN,fin,NaN,NaN,0,"nufft");









    



Single pseudo Wigner Ville
Use nufft.



In [12]:

    
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.tfrshow(abs.(tfrpf),xs[2]-xs[1],xs[1],xs[end],NaN,NaN,0.7*130,"Spectral")
PyPlot.xlabel("time [day]")
PyPlot.title("Zoom-up (FFT)")
PyPlot.ylim(0.98,1.02)
ax = fig[:add_subplot](1,2,2)
a=juwplot.tfrshow(abs.(tfrpfn),xs[2]-xs[1],xs[1],xs[end],fin[1],fin[end],0.7,"Spectral")
PyPlot.xlabel("time [day]")
PyPlot.title("NuFFT")









    












    Out[12]:





PyObject Text(0.5,1,'NuFFT')



In [13]:

    
indfn=extif.maxif(abs.(tfrpfn));
fn=juwutils.index_to_frequency(indfn, fin,xs[2]-xs[1],nsample,nnufft);

comparison of the extracted IFs

the red curve corresponds to the input IF, gray is the pseudo Wigner-Ville with FFT, and blue is the pseudo Wigner-Ville with NUFFT.
The NUFFT is necessary.
The pseudo Wigner-Ville with NUFFT provides us the IF with the sufficient resolution.



In [14]:

    
PyPlot.title("no noise")
PyPlot.xlabel("time")
PyPlot.ylabel("instantaneous frequency")
PyPlot.ylim(0.995,1.005)
PyPlot.plot(xs,fx, color="gray",lw=3)
PyPlot.plot(xs,fn, color="C0",lw=3)
PyPlot.plot(xs,iws/(2*pi), color="C3", ls="dashed",lw=3)
PyPlot.legend(["IF from PWV (FFT)","IF from PWV (NUFFT)","Input IF"])









    












    Out[14]:





PyObject <matplotlib.legend.Legend object at 0x7fa8d8933ba8>

Thinning out time grids for reducing computational time



In [15]:

    
# using one in 20 time grids 
nthin=20
itc=collect(1:nthin:nsample);
tfrpfi=cohenclass.tfrpwv(z,NaN,NaN,fin,itc,NaN,0,"nufft");









    



Single pseudo Wigner Ville
Use nufft.



In [16]:

    
fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.tfrshow(abs.(tfrpfn),xs[2]-xs[1],xs[1],xs[end],fin[1],fin[end],0.7,"Spectral")
PyPlot.xlabel("time [day]")
PyPlot.title("NuFFT")
ax = fig[:add_subplot](1,2,2)
a=juwplot.tfrshow(abs.(tfrpfi),(xs[2]-xs[1])*nthin,xs[1],xs[end],fin[1],fin[end],0.7,"Spectral")
PyPlot.xlabel("time [day]")
PyPlot.title("NuFFT & thinning")









    












    Out[16]:





PyObject Text(0.5,1,'NuFFT & thinning')

this technique is useful for large dataset



In [ ]: