19
19
assumption is that the sampling rate is 2kHz.
20
20
21
21
"""
22
- function cwt (Y:: AbstractArray{T,N} , cWav:: CWT , daughters, fftPlans = 1 ) where {N, T}
23
- @assert typeof (N)<: Integer
22
+ function cwt (Y:: AbstractArray{T,N} , cWav:: CWT , daughters, fftPlans = 1 ) where {N,T}
23
+ @assert typeof (N) <: Integer
24
24
# vectors behave a bit strangely, so we reshape them
25
- if N== 1
26
- Y= reshape (Y,(length (Y), 1 ))
25
+ if N == 1
26
+ Y = reshape (Y, (length (Y), 1 ))
27
27
end
28
28
n1 = size (Y, 1 )
29
29
@@ -35,8 +35,8 @@ function cwt(Y::AbstractArray{T,N}, cWav::CWT, daughters, fftPlans = 1) where {N
35
35
x̂, fromPlan = prepSignalAndPlans (x, cWav, fftPlans)
36
36
# If the vector isn't long enough to actually have any other scales, just
37
37
# return the averaging
38
- if nScales <= 0 || size (daughters,2 ) == 1
39
- daughters = daughters[:,1 : 1 ]
38
+ if nScales <= 0 || size (daughters, 2 ) == 1
39
+ daughters = daughters[:, 1 : 1 ]
40
40
nScales = 1
41
41
end
42
42
@@ -61,8 +61,8 @@ function cwt(Y::AbstractArray{T,N}, cWav::CWT, daughters, fftPlans = 1) where {N
61
61
wave = permutedims (wave, [1 , ndims (wave), (2 : (ndims (wave)- 1 )). .. ])
62
62
ax = axes (wave)
63
63
wave = wave[1 : n1, ax[2 : end ]. .. ]
64
- if N== 1
65
- wave = dropdims (wave, dims= 3 )
64
+ if N == 1
65
+ wave = dropdims (wave, dims = 3 )
66
66
end
67
67
68
68
return wave
82
82
# yes | both | fft
83
83
# no | rfft | fft
84
84
# Analytic on Real input
85
- function prepSignalAndPlans (x:: AbstractArray{T} , cWav:: CWT{W,S,WaTy, N, true} , fftPlans) where {T <: Real , W,S,WaTy,N}
85
+ function prepSignalAndPlans (x:: AbstractArray{T} , cWav:: CWT{W,S,WaTy,N, true} , fftPlans) where {T<: Real ,W,S,WaTy,N}
86
86
# analytic wavelets that are being applied on real inputs
87
- if fftPlans isa Tuple{<: AbstractFFTs.Plan{<:Real} , <: AbstractFFTs.Plan{<:Complex} }
87
+ if fftPlans isa Tuple{<: AbstractFFTs.Plan{<:Real} ,<: AbstractFFTs.Plan{<:Complex} }
88
88
# they handed us the right kind of thing, so no need to make new ones
89
89
x̂ = fftPlans[1 ] * x
90
90
fromPlan = fftPlans[2 ]
91
91
else
92
- toPlan = plan_rfft (x,1 )
92
+ toPlan = plan_rfft (x, 1 )
93
93
x̂ = toPlan * x
94
- fromPlan = plan_fft (x,1 )
94
+ fromPlan = plan_fft (x, 1 )
95
95
end
96
96
return x̂, fromPlan
97
97
end
98
98
99
99
# Non-analytic on Real input
100
- function prepSignalAndPlans (x:: AbstractArray{T} , cWav:: CWT{W,S,WaTy, N, false} , fftPlans) where {T <: Real , W,S,WaTy,N}
100
+ function prepSignalAndPlans (x:: AbstractArray{T} , cWav:: CWT{W,S,WaTy,N, false} , fftPlans) where {T<: Real ,W,S,WaTy,N}
101
101
# real wavelets that are being applied on real inputs
102
102
if fftPlans isa AbstractFFTs. Plan{<: Real }
103
103
# they handed us the right kind of thing, so no need to make new ones
104
104
x̂ = fftPlans * x
105
105
fromPlan = fftPlans
106
106
else
107
- fromPlan = plan_rfft (x,1 )
107
+ fromPlan = plan_rfft (x, 1 )
108
108
x̂ = fromPlan * x
109
109
end
110
110
return x̂, fromPlan
111
111
end
112
112
# complex input
113
- function prepSignalAndPlans (x:: AbstractArray{T} , cWav, fftPlans) where {T <: Complex , W,S,WaTy,N }
113
+ function prepSignalAndPlans (x:: AbstractArray{T} , cWav, fftPlans) where {T<: Complex }
114
114
# real wavelets that are being applied on real inputs
115
115
if fftPlans isa AbstractFFTs. Plan{<: Complex }
116
116
# they handed us the right kind of thing, so no need to make new ones
117
117
x̂ = fftPlans * x
118
118
fromPlan = fftPlans
119
119
else
120
- fromPlan = plan_fft (x,1 )
120
+ fromPlan = plan_fft (x, 1 )
121
121
x̂ = fromPlan * x
122
122
end
123
123
return x̂, fromPlan
124
124
end
125
125
126
126
# analytic on real data with an averaging function
127
- function analyticTransformReal! (wave, daughters, x̂, fftPlan, averagingType :: Union{Father, Dirac} )
127
+ function analyticTransformReal! (wave, daughters, x̂, fftPlan, :: Union{Father,Dirac} )
128
128
outer = axes (x̂)[2 : end ]
129
129
n1 = size (x̂, 1 )
130
- isSourceEven = mod (size (wave,1 ) + 1 , 2 )
130
+ isSourceEven = mod (size (wave, 1 ) + 1 , 2 )
131
131
# the averaging function isn't analytic, so we need to do both positive and
132
132
# negative frequencies
133
- @views tmpWave = x̂ .* daughters[:,1 ]
134
- @views wave[(n1+ 1 ): end , outer... , 1 ] = reverse (conj .(tmpWave[2 : end - isSourceEven, outer... ]), dims= 1 )
133
+ @views tmpWave = x̂ .* daughters[:, 1 ]
134
+ @views wave[(n1+ 1 ): end , outer... , 1 ] = reverse (conj .(tmpWave[2 : end - isSourceEven, outer... ]), dims = 1 )
135
135
@views wave[1 : n1, outer... , 1 ] = tmpWave
136
136
@views wave[:, outer... , 1 ] = fftPlan \ (wave[:, outer... , 1 ]) # averaging
137
- for j in 2 : size (daughters,2 )
138
- @views wave[1 : n1, outer... , j] = x̂ .* daughters[:,j]
137
+ for j = 2 : size (daughters, 2 )
138
+ @views wave[1 : n1, outer... , j] = x̂ .* daughters[:, j]
139
139
wave[:, outer... , j] = fftPlan \ (wave[:, outer... , j]) # wavelet transform
140
140
end
141
141
end
142
142
143
143
# analytic on complex data with an averaging function
144
- function analyticTransformComplex! (wave, daughters, x̂, fftPlan, averagingType :: Union{Father, Dirac} )
144
+ function analyticTransformComplex! (wave, daughters, x̂, fftPlan, :: Union{Father,Dirac} )
145
145
outer = axes (x̂)[2 : end ]
146
146
n1 = size (daughters, 1 )
147
- isSourceEven = mod (size (wave,1 ) + 1 , 2 )
147
+ isSourceEven = mod (size (wave, 1 ) + 1 , 2 )
148
148
# the averaging function isn't analytic, so we need to do both positive and
149
149
# negative frequencies
150
- @views positiveFreqs = x̂[1 : n1, outer... ] .* daughters[:,1 ]
151
- @views negativeFreqs = x̂[(n1- isSourceEven+ 1 ): end , outer... ] .* reverse (conj .(daughters[2 : end ,1 ]))
150
+ @views positiveFreqs = x̂[1 : n1, outer... ] .* daughters[:, 1 ]
151
+ @views negativeFreqs = x̂[(n1- isSourceEven+ 1 ): end , outer... ] .* reverse (conj .(daughters[2 : end , 1 ]))
152
152
@views wave[(n1- isSourceEven+ 1 ): end , outer... , 1 ] = negativeFreqs
153
153
@views wave[1 : n1, outer... , 1 ] = positiveFreqs
154
154
@views wave[:, outer... , 1 ] = fftPlan \ (wave[:, outer... , 1 ]) # averaging
155
- for j in 2 : size (daughters,2 )
156
- @views wave[1 : n1, outer... , j] = x̂[1 : n1, outer... ] .* daughters[:,j]
155
+ for j = 2 : size (daughters, 2 )
156
+ @views wave[1 : n1, outer... , j] = x̂[1 : n1, outer... ] .* daughters[:, j]
157
157
@views wave[:, outer... , j] = fftPlan \ (wave[:, outer... , j]) # wavelet transform
158
158
end
159
159
end
160
160
161
161
function analyticTransformComplex! (wave, daughters, x̂, fftPlan, averagingType)
162
162
outer = axes (x̂)[2 : end ]
163
163
n1 = size (x̂, 1 )
164
- for j in 1 : size (daughters,2 )
165
- @views wave[1 : n1, outer... , j] = x̂[1 : n1, outer... ] .* daughters[:,j]
164
+ for j = 1 : size (daughters, 2 )
165
+ @views wave[1 : n1, outer... , j] = x̂[1 : n1, outer... ] .* daughters[:, j]
166
166
@views wave[:, outer... , j] = fftPlan \ (wave[:, outer... , j]) # wavelet transform
167
167
end
168
168
end
169
169
170
170
# analytic on real data without an averaging function
171
- function analyticTransformReal! (wave, daughters, x̂, fftPlan, averagingType :: NoAve )
171
+ function analyticTransformReal! (wave, daughters, x̂, fftPlan, :: NoAve )
172
172
outer = axes (x̂)[2 : end ]
173
173
n1 = size (x̂, 1 )
174
174
# the no averaging version
175
- for j in 1 : size (daughters,2 )
176
- wave[1 : n1, outer... , j] = x̂ .* daughters[:,j]
175
+ for j = 1 : size (daughters, 2 )
176
+ wave[1 : n1, outer... , j] = x̂ .* daughters[:, j]
177
177
wave[:, outer... , j] = fftPlan \ (wave[:, outer... , j]) # wavelet transform
178
178
end
179
179
end
@@ -183,8 +183,8 @@ function otherwiseTransform!(wave::AbstractArray{<:Real}, daughters, x̂, fromPl
183
183
# real wavelets on real data: that just makes sense
184
184
outer = axes (x̂)[2 : end ]
185
185
n1 = size (x̂, 1 )
186
- for j in 1 : size (daughters, 2 )
187
- @views tmp = x̂ .* daughters[:,j]
186
+ for j = 1 : size (daughters, 2 )
187
+ @views tmp = x̂ .* daughters[:, j]
188
188
@views wave[:, outer... , j] = fromPlan \ tmp # wavelet transform
189
189
end
190
190
end
@@ -195,35 +195,34 @@ function otherwiseTransform!(wave::AbstractArray{<:Complex}, daughters, x̂, fro
195
195
outer = axes (x̂)[2 : end ]
196
196
n1 = size (daughters, 1 )
197
197
isSourceEven = mod (size (fromPlan, 1 ) + 1 , 2 )
198
- for j in 1 : size (daughters, 2 )
199
- @views wave[1 : n1, outer... , j] = @views x̂[1 : n1,outer... ] .* daughters[:,j]
200
- @views wave[n1- isSourceEven+ 1 : end , outer... , j] = x̂[n1- isSourceEven+ 1 : end ,outer... ] .* reverse (conj .(daughters[2 : end ,j]))
198
+ for j = 1 : size (daughters, 2 )
199
+ @views wave[1 : n1, outer... , j] = @views x̂[1 : n1, outer... ] .* daughters[:, j]
200
+ @views wave[n1- isSourceEven+ 1 : end , outer... , j] = x̂[n1- isSourceEven+ 1 : end , outer... ] .* reverse (conj .(daughters[2 : end , j]))
201
201
@views wave[:, outer... , j] = fromPlan \ (wave[:, outer... , j]) # wavelet transform
202
202
end
203
203
end
204
204
205
205
function reflect (Y, bt)
206
206
n1 = size (Y, 1 )
207
207
if typeof (bt) <: ZPBoundary
208
- base2 = ceil (Int, log2 (n1)); # power of 2 nearest to N
209
- x = cat (Y, zeros (2 ^ (base2)- n1, size (Y)[2 : end ]. .. ), dims= 1 )
208
+ base2 = ceil (Int, log2 (n1)) # power of 2 nearest to N
209
+ x = cat (Y, zeros (2 ^ (base2) - n1, size (Y)[2 : end ]. .. ), dims = 1 )
210
210
elseif typeof (bt) <: SymBoundary
211
- x = cat (Y, reverse (Y,dims = 1 ), dims = 1 )
211
+ x = cat (Y, reverse (Y, dims = 1 ), dims = 1 )
212
212
else
213
213
x = Y
214
214
end
215
215
return x
216
216
end
217
217
218
218
219
- function cwt (Y:: AbstractArray{T} , c:: CWT{W} ; varArgs... ) where {T<: Number , S<: Real , V<: Real ,
220
- W<: WaveletBoundary }
221
- daughters,ω = computeWavelets (size (Y, 1 ), c; varArgs... )
219
+ function cwt (Y:: AbstractArray{T} , c:: CWT{W} ; varArgs... ) where {T<: Number ,W<: WaveletBoundary }
220
+ daughters, ω = computeWavelets (size (Y, 1 ), c; varArgs... )
222
221
return cwt (Y, c, daughters)
223
222
end
224
223
225
- cwt (Y:: AbstractArray{T} , w:: ContWaveClass ; varargs... ) where {T<: Number , S <: Real , V <: Real } = cwt (Y,CWT (w); varargs... )
226
- cwt (Y:: AbstractArray{T} ) where T<: Real = cwt (Y,Morlet ())
224
+ cwt (Y:: AbstractArray{T} , w:: ContWaveClass ; varargs... ) where {T<: Number } = cwt (Y, CWT (w); varargs... )
225
+ cwt (Y:: AbstractArray{T} ) where { T<: Real } = cwt (Y, Morlet ())
227
226
228
227
abstract type InverseType end
229
228
struct DualFrames <: InverseType end
@@ -244,35 +243,35 @@ Return the inverse continuous wavelet transform, computed using the simple dual
244
243
Return the inverse continuous wavelet transform, computed using the canonical dual frame ``\\ tilde{\\ widehat{ψ}} = \\ frac{ψ̂_n(ω)}{∑_n\\ |ψ̂_n(ω)\\ |^2}``. The algorithm is to compute the cwt again, but using the canonical dual frame; consiquentially, it is the most computationally intensive of the three algorithms, and typically the best behaved. Will be numerically unstable if the high frequencies of all of the wavelets are too small however, and tends to fail spectacularly in this case.
245
244
246
245
"""
247
- function icwt (res:: AbstractArray , cWav:: CWT , inverseStyle :: PenroseDelta )
248
- Ŵ = computeWavelets (size (res,1 ), cWav)[1 ]
246
+ function icwt (res:: AbstractArray , cWav:: CWT , :: PenroseDelta )
247
+ Ŵ = computeWavelets (size (res, 1 ), cWav)[1 ]
249
248
β = computeDualWeights (Ŵ, cWav)
250
249
testDualCoverage (β, Ŵ)
251
- compXRecon = sum (res .* β, dims= 2 )
252
- imagXRecon = irfft (im* rfft (imag .(compXRecon),1 ), size (compXRecon,1 )) # turns out the dual frame for the imaginary part is rather gross in the time domain
250
+ compXRecon = sum (res .* β, dims = 2 )
251
+ imagXRecon = irfft (im * rfft (imag .(compXRecon), 1 ), size (compXRecon, 1 )) # turns out the dual frame for the imaginary part is rather gross in the time domain
253
252
return imagXRecon + real .(compXRecon)
254
253
end
255
254
256
- function icwt (res:: AbstractArray , cWav:: CWT , inverseStyle :: NaiveDelta )
257
- Ŵ = computeWavelets (size (res,1 ), cWav)[1 ]
258
- β = computeNaiveDualWeights (Ŵ, cWav, size (res,1 ))
255
+ function icwt (res:: AbstractArray , cWav:: CWT , :: NaiveDelta )
256
+ Ŵ = computeWavelets (size (res, 1 ), cWav)[1 ]
257
+ β = computeNaiveDualWeights (Ŵ, cWav, size (res, 1 ))
259
258
testDualCoverage (β, Ŵ)
260
- compXRecon = sum (res .* β, dims= 2 )
261
- imagXRecon = irfft (im* rfft (imag .(compXRecon),1 ), size (compXRecon,1 )) # turns out the dual frame for the imaginary part is rather gross in the time domain
259
+ compXRecon = sum (res .* β, dims = 2 )
260
+ imagXRecon = irfft (im * rfft (imag .(compXRecon), 1 ), size (compXRecon, 1 )) # turns out the dual frame for the imaginary part is rather gross in the time domain
262
261
return imagXRecon + real .(compXRecon)
263
262
end
264
263
265
- function icwt (res:: AbstractArray , cWav:: CWT , inverseStyle :: DualFrames )
266
- Ŵ = computeWavelets (size (res,1 ), cWav)[1 ]
264
+ function icwt (res:: AbstractArray , cWav:: CWT , :: DualFrames )
265
+ Ŵ = computeWavelets (size (res, 1 ), cWav)[1 ]
267
266
canonDualFrames = explicitConstruction (Ŵ)
268
267
testDualCoverage (canonDualFrames)
269
268
convolved = cwt (res, cWav, canonDualFrames)
270
269
ax = axes (convolved)
271
- @views xRecon = sum (convolved[:,i,i, ax[4 : end ]. .. ] for i= 1 : size (Ŵ,2 ))
270
+ @views xRecon = sum (convolved[:, i, i, ax[4 : end ]. .. ] for i = 1 : size (Ŵ, 2 ))
272
271
return xRecon
273
272
end
274
273
275
- icwt (Y:: AbstractArray , w:: ContWaveClass ; varargs... ) where {S <: Real , T <: Real , V <: Real } = icwt (Y,CWT (w), PenroseDelta (); varargs... )
276
- icwt (Y:: AbstractArray ; varargs... ) = icwt (Y,Morlet (), PenroseDelta (); varargs... )
274
+ icwt (Y:: AbstractArray , w:: ContWaveClass ; varargs... ) = icwt (Y, CWT (w), PenroseDelta (); varargs... )
275
+ icwt (Y:: AbstractArray ; varargs... ) = icwt (Y, Morlet (), PenroseDelta (); varargs... )
277
276
278
277
# CWT (continuous wavelet transform directly) TODO : direct if sufficiently small
0 commit comments