Skip to content

Commit

Permalink
Merge df2d2d7 into 140cb53
Browse files Browse the repository at this point in the history
  • Loading branch information
mtaoraki committed Oct 7, 2015
2 parents 140cb53 + df2d2d7 commit 861d000
Show file tree
Hide file tree
Showing 3 changed files with 144 additions and 46 deletions.
19 changes: 19 additions & 0 deletions example/continuoustransform1d.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
using Wavelets
using PyPlot

dt=0.01;
t=1e-5:dt:40;
Y=0.1*randn(length(t))+sin(2*π*5*t);

wave,period,scale,coi=Wavelets.cwtft(Y,dt);

f=1./period;
ax=subplot(111);
cf=ax[:contourf](t,f,log(float(abs(wave).^2)),50);
cb=colorbar(cf);
cb[:set_label]("Wavelet Power");
ax[:fill_between](t,1./coi, color="white", alpha=0.5);
yscale("log");
xlabel("Time [s]");
ylabel("Frequency [Hz]");
ylim([minimum(f),maximum(f)]);
112 changes: 77 additions & 35 deletions src/transforms.jl
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
module Transforms
using ..Util, ..WT
using Compat
export dwt, idwt, dwt!, idwt!,
wpt, iwpt, wpt!, iwpt!,
dwtc, idwtc
export dwt, idwt, dwt!, idwt!,
wpt, iwpt, wpt!, iwpt!,
dwtc, idwtc, cwtft
VERSION < v"0.4-" && using Docile

# TODO Use StridedArray instead of AbstractArray where writing to array.
Expand All @@ -13,31 +13,31 @@ typealias WPTArray AbstractVector

# DWT (discrete wavelet transform)

for (Xwt, Xwt!, _Xwt!, fw) in ((:dwt, :dwt!, :_dwt!, true),
for (Xwt, Xwt!, _Xwt!, fw) in ((:dwt, :dwt!, :_dwt!, true),
(:idwt, :idwt!, :_dwt!, false))
@eval begin
# filter
function ($Xwt){T<:FloatingPoint}(x::DWTArray{T},
filter::OrthoFilter,
function ($Xwt){T<:FloatingPoint}(x::DWTArray{T},
filter::OrthoFilter,
L::Integer=maxtransformlevels(x))
y = Array(T, size(x))
return ($_Xwt!)(y, x, filter, L, $fw)
end
function ($Xwt!){T<:FloatingPoint}(y::DWTArray{T}, x::DWTArray{T},
filter::OrthoFilter,
function ($Xwt!){T<:FloatingPoint}(y::DWTArray{T}, x::DWTArray{T},
filter::OrthoFilter,
L::Integer=maxtransformlevels(x))
return ($_Xwt!)(y, x, filter, L, $fw)
end
# lifting
function ($Xwt){T<:FloatingPoint}(x::DWTArray{T},
scheme::GLS,
function ($Xwt){T<:FloatingPoint}(x::DWTArray{T},
scheme::GLS,
L::Integer=maxtransformlevels(x))
y = Array(T, size(x))
copy!(y, x)
return ($_Xwt!)(y, scheme, L, $fw)
end
function ($Xwt!){T<:FloatingPoint}(y::DWTArray{T},
scheme::GLS,
function ($Xwt!){T<:FloatingPoint}(y::DWTArray{T},
scheme::GLS,
L::Integer=maxtransformlevels(x))
return ($_Xwt!)(y, scheme, L, $fw)
end
Expand Down Expand Up @@ -99,46 +99,46 @@ The inverse of `dwt!`.

# WPT (wavelet packet transform)

for (Xwt, Xwt!, _Xwt!, fw) in ((:wpt, :wpt!, :_wpt!, true),
for (Xwt, Xwt!, _Xwt!, fw) in ((:wpt, :wpt!, :_wpt!, true),
(:iwpt, :iwpt!, :_wpt!, false))
@eval begin
function ($Xwt){T<:FloatingPoint}(x::WPTArray{T},
wt::DiscreteWavelet,
function ($Xwt){T<:FloatingPoint}(x::WPTArray{T},
wt::DiscreteWavelet,
L::Integer=maxtransformlevels(x))
return ($Xwt)(x, wt, maketree(length(x), L, :full))
end
# filter
function ($Xwt){T<:FloatingPoint}(x::WPTArray{T},
filter::OrthoFilter,
function ($Xwt){T<:FloatingPoint}(x::WPTArray{T},
filter::OrthoFilter,
tree::BitVector=maketree(x, :full))
y = Array(T, size(x))
return ($_Xwt!)(y, x, filter, tree, $fw)
end
function ($Xwt!){T<:FloatingPoint}(y::WPTArray{T}, x::WPTArray{T},
filter::OrthoFilter,
function ($Xwt!){T<:FloatingPoint}(y::WPTArray{T}, x::WPTArray{T},
filter::OrthoFilter,
tree::BitVector=maketree(x, :full))
return ($_Xwt!)(y, x, filter, tree, $fw)
end
function ($Xwt!){T<:FloatingPoint}(y::WPTArray{T}, x::WPTArray{T},
filter::OrthoFilter,
function ($Xwt!){T<:FloatingPoint}(y::WPTArray{T}, x::WPTArray{T},
filter::OrthoFilter,
L::Integer=maxtransformlevels(x))
return ($Xwt!)(y, x, filter, maketree(length(x), L, :full))
end
# lifting
function ($Xwt){T<:FloatingPoint}(x::WPTArray{T},
scheme::GLS,
function ($Xwt){T<:FloatingPoint}(x::WPTArray{T},
scheme::GLS,
tree::BitVector=maketree(x, :full))
y = Array(T, size(x))
copy!(y, x)
return ($_Xwt!)(y, scheme, tree, $fw)
end
function ($Xwt!){T<:FloatingPoint}(y::WPTArray{T},
scheme::GLS,
function ($Xwt!){T<:FloatingPoint}(y::WPTArray{T},
scheme::GLS,
tree::BitVector=maketree(y, :full))
return ($_Xwt!)(y, scheme, tree, $fw)
end
function ($Xwt!){T<:FloatingPoint}(y::WPTArray{T},
scheme::GLS,
function ($Xwt!){T<:FloatingPoint}(y::WPTArray{T},
scheme::GLS,
L::Integer=maxtransformlevels(x))
return ($Xwt!)(y, scheme, maketree(length(x), L, :full))
end
Expand All @@ -159,7 +159,50 @@ end # for
#icwt(::AbstractVector, ::ContinuousWavelet)

# CWTFT (continuous wavelet transform via FFT)
#cwtft(::AbstractVector, ::ContinuousWavelet)

## After Torrence & Compo (1998):

function cwtft{T<:Real}(Y::AbstractArray{T},dt::Number; pad::Bool=false,dj::Number=0.25,s0::Number=2.*dt,J1::Number=-1,mother::WT.ContinuousWavelet=WT.morlet,param::Number=WT.sparam(mother))

#Y=Y[:];
n1 = length(Y);

if J1 == -1
J1=floor(Int,(log(n1*dt/s0)/log(2.))/dj);
end
#....construct time series to analyze, pad if necessary
x = Y - mean(Y);
if pad
base2 = floor(Int,log(n1)/log(2) + 0.4999); # power of 2 nearest to N
x = [x, zeros(2^(floor(Int,base2)+1)-n1)];
end
n = length(x);

#....construct wavenumber array used in transform [Eqn(5)]
k = 1:floor(Int,n/2);
k = [0.; k; -k[floor(Int,(n-1)/2):-1:1]]*((2*pi)/(n*dt));
#....compute FFT of the (padded) time series
f = fft(x); # [Eqn(3)]
#....construct SCALE array & empty PERIOD & WAVE arrays
scale = s0*2.^((0:J1)*dj);

wave = zeros(Complex,J1+1,n); # define the wavelet array


# loop through all scales and compute transform
for a1 in 1:J1+1
daughter=WT.Daughter(mother,scale[a1],k,param,n)
wave[a1,:] = ifft(f.*daughter) # wavelet transform[Eqn(4)]
end
fourier_factor=WT.FourierFactor(mother,param);
period = fourier_factor*scale;
coi = WT.COI(mother,fourier_factor).*dt*[1E-5; 1:((n1+1)/2-1); flipdim((1:(n1/2-1)),1); 1E-5]; # COI [Sec.3g]
wave = wave[:,1:n1]; # get rid of padding before returning


return wave,period,scale,coi
end

#icwtft(::AbstractVector, ::ContinuousWavelet)

@deprecate dwt!(y, x, filter::OrthoFilter, L, fw) (fw ? dwt!(y,x,filter,L) : idwt!(y,x,filter,L))
Expand All @@ -178,14 +221,14 @@ end
for (Xwt_oop!, Xwt!) in ((:dwt_oop!, :dwt!), (:idwt_oop!, :idwt!))
@eval begin
# filter
function ($Xwt_oop!){T<:FloatingPoint}(y::DWTArray{T}, x::DWTArray{T},
filter::OrthoFilter,
function ($Xwt_oop!){T<:FloatingPoint}(y::DWTArray{T}, x::DWTArray{T},
filter::OrthoFilter,
L::Integer=maxtransformlevels(x))
return ($Xwt!)(y, x, filter, L)
end
# lifting
function ($Xwt_oop!){T<:FloatingPoint}(y::DWTArray{T}, x::DWTArray{T},
scheme::GLS,
function ($Xwt_oop!){T<:FloatingPoint}(y::DWTArray{T}, x::DWTArray{T},
scheme::GLS,
L::Integer=maxtransformlevels(x))
copy!(y, x)
return ($Xwt!)(y, scheme, L)
Expand All @@ -207,12 +250,12 @@ for (Xwtc, Xwt) in ((:dwtc, :dwt!), (:idwtc, :idwt!))
y = Array(T, sizex)
xc = Array(T, sizexc)
yc = Array(T, sizexc)

ind = Array(Any, dim)
for i = 1:dim
ind[i] = 1:sizex[i]
end

for d = 1:sizex[td]
ind[td] = d
if dim > 2
Expand Down Expand Up @@ -267,4 +310,3 @@ include("transforms_lifting.jl")


end # module

59 changes: 48 additions & 11 deletions src/wt.jl
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
module WT
export DiscreteWavelet,
ContinuousWavelet,
FilterWavelet,
LSWavelet,
OrthoFilter,
Expand All @@ -14,7 +15,7 @@ VERSION < v"0.4-" && using Docile
# TYPE HIERARCHY

abstract DiscreteWavelet{T}
#abstract ContinuousWavelet{T}
abstract ContinuousWavelet
# discrete transforms via filtering
abstract FilterWavelet{T} <: DiscreteWavelet{T}
# discrete transforms via lifting
Expand Down Expand Up @@ -128,7 +129,7 @@ end
Wavelet type for discrete orthogonal transforms by filtering.
**See also:** `GLS`, `wavelet`
""" ->
""" ->
immutable OrthoFilter{T<:WaveletBoundary} <: FilterWavelet{T}
qmf ::Vector{Float64} # quadrature mirror filter
name ::ASCIIString # filter short name
Expand Down Expand Up @@ -205,11 +206,11 @@ Base.length(s::LSStep) = length(s.param)
Base.length(s::LSStepParam) = length(s.coef)

@doc """
Wavelet type for discrete general (bi)orthogonal transforms
Wavelet type for discrete general (bi)orthogonal transforms
by using a lifting scheme.
**See also:** `OrthoFilter`, `wavelet`
""" ->
""" ->
immutable GLS{T<:WaveletBoundary} <: LSWavelet{T}
step ::Vector{LSStep{Float64}} # steps to be taken
norm1 ::Float64 # normalization of scaling coefs.
Expand All @@ -229,7 +230,45 @@ name(s::GLS) = s.name

# IMPLEMENTATIONS OF ContinuousWavelet

# ...
for (TYPE, NAMEBASE, PARAM, FFFUNCTION, COIFUNCTION, DFUNCTION, SUPERCLASS) in (
(:Morlet, "morlet",
6., #PARAM
:((4*π)/(m + sqrt(2 + m^2))), #FourierFactorFunction
:(1/sqrt(2)), #COIFunction
:(daughter=sqrt(scale*k[2])*(π^(-1//4))*sqrt(n)*exp(-(scale*k - m).^2/2)), #DaughterFunction
:ContinuousWavelet),
(:Paul, "paul",
4., #PARAM
:(4*pi/(2*m+1)), #FourierFactorFunction
:(1*sqrt(2)), #COIFunction
:(daughter=sqrt(scale*k[2])*(2^m/sqrt(m*prod(2:(2*m-1))))*sqrt(n)*exp(-(scale*k))), #DaughterFunction
:ContinuousWavelet),
(:DOG, "dog",
2., #PARAM
:(2*pi*sqrt(2./(2*m+1))), #FourierFactorFunction
:(1/sqrt(2)), #COIFunction
:(expnt = -(scale*k).^2/2.0;
norm = sqrt(scale*k[2]/gamma(m+0.5))*sqrt(n);
daughter = -norm*(i^m)*((scale*k2)^m).*exp(expnt)), #DaughterFunction
:ContinuousWavelet) # TODO moments
)
@eval begin
immutable $TYPE <: $SUPERCLASS end
name(::$TYPE) = string($NAMEBASE)::ASCIIString
sparam(::$TYPE) = $PARAM
FourierFactor(::$TYPE,m::Real)= $FFFUNCTION
COI(::$TYPE,FFactor::Real)=FFactor* $COIFUNCTION
function Daughter(::$TYPE,scale::Real,k::Array{Float64,1},m::Real,n::Integer)
$DFUNCTION
daughter[k.<0]=0.;
return daughter
end
end
CONSTNAME = symbol(NAMEBASE)
@eval begin
const $CONSTNAME = $TYPE() # type shortcut
end
end

# TRANSFORM TYPE CONSTRUCTORS

Expand Down Expand Up @@ -278,7 +317,7 @@ function daubechies(N::Int)
# Find roots in y domain (truncated binomial series; (1 - y)^{-N})
Y = roots(C)

# Find roots in z domain:
# Find roots in z domain:
# z + z^{-1} = 2 - 4*y
# where y is a root from above
Z = zeros(Complex128, 2*N-2)
Expand Down Expand Up @@ -347,7 +386,7 @@ function vieta(R::AbstractVector)
C[1] = 1
Ci::Complex128 = 0
Cig::Complex128 = 0

@inbounds for k = 1:n
Ci = C[1]
for i = 1:k
Expand All @@ -361,7 +400,7 @@ end


# scaling filters h (low pass)
# the number at end of a filter name is the
# the number at end of a filter name is the
# number of vanishing moments of the mother wavelet function
# sources:
# http://statweb.stanford.edu/~wavelab/ (Orthogonal/MakeONFilter.m)
Expand Down Expand Up @@ -476,11 +515,9 @@ const SCHEMES = @compat Dict{ASCIIString,NTuple{3}}(
0.5176380902050414,
1.9318516525781364)

)
)




end


0 comments on commit 861d000

Please sign in to comment.