Skip to content

Commit c11fe30

Browse files
bjarthurbjarthur
committed
add F-test to mt_pgram
1 parent 08f4b0f commit c11fe30

File tree

4 files changed

+71
-17
lines changed

4 files changed

+71
-17
lines changed

‎Project.toml‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@ version = "0.8.3"
44

55
[deps]
66
Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
7+
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
78
FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
89
IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
910
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -22,6 +23,7 @@ OffsetArraysExt = "OffsetArrays"
2223
[compat]
2324
Bessels = "0.2"
2425
DelimitedFiles = "1.6"
26+
Distributions = "0.25.119"
2527
FFTW = "1.8"
2628
IterTools = "1.4"
2729
LinearAlgebra = "1.6"

‎src/multitaper.jl‎

Lines changed: 39 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,8 @@
22
##### Multitapered periodogram
33
#####
44

5+
using Distributions
6+
57
struct MTConfig{T,R1,F,P,T1,T2,W,R2}
68
n_samples::Int
79
fs::R1
@@ -13,11 +15,12 @@ struct MTConfig{T,R1,F,P,T1,T2,W,R2}
1315
fft_output_tmp::T2
1416
window::W
1517
onesided::Bool
18+
ftest::Bool
1619
r::R2 # inverse normalization; e.g. equal to `fs*N` for an unwindowed/untapered periodogram of a signal of length `N`
1720
# e.g. equal to `fs*ntapers*ones(ntapers)` when tapered by properly normalized window functions (i.e. norm-2 equal to 1 for each taper)
1821
# can be adjusted to weight the tapers, e.g. by the eigenvalues of the DPSS windows
1922
function MTConfig{T}(n_samples, nfft, ntapers, freq::F, fs::R1, plan::P, fft_input_tmp::T1,
20-
fft_output_tmp::T2, window::W, onesided, r::R2) where {T,R1,F,P,T1,T2,W,R2}
23+
fft_output_tmp::T2, window::W, onesided, ftest, r::R2) where {T,R1,F,P,T1,T2,W,R2}
2124
n_samples > 0 || throw(ArgumentError("`n_samples` must be positive"))
2225
nfft >= n_samples || throw(ArgumentError("Must have `nfft >= n_samples`"))
2326
ntapers > 0 || throw(ArgumentError("`ntapers` must be positive"))
@@ -44,7 +47,7 @@ struct MTConfig{T,R1,F,P,T1,T2,W,R2}
4447
ntapers, freq, plan,
4548
fft_input_tmp,
4649
fft_output_tmp, window,
47-
onesided, r)
50+
onesided, ftest, r)
4851
end
4952
end
5053

@@ -85,6 +88,7 @@ end
8588
ntapers = 2 * nw - 1,
8689
taper_weights = fill(1/ntapers, ntapers),
8790
onesided::Bool=T<:Real,
91+
ftest::Bool=false,
8892
fft_flags = FFTW.MEASURE)
8993
9094
Creates a config object which holds the configuration state
@@ -107,11 +111,12 @@ as none of the input arguments change.
107111
The default setting is to simply average them.
108112
* `onesided`: whether or not to compute a "one-sided" FFT by using that real signal data
109113
yields conjugate-symmetry in Fourier space.
114+
* `ftest`: whether or not to compute the significance of harmonic line compenents as described in Thomson (1982; 10.1109/PROC.1982.12433)
110115
* `fft_flags`: flags to control how the FFT plan is generated.
111116
"""
112117
function MTConfig{T}(n_samples; fs=1, nfft=nextpow(2, n_samples), window=nothing, nw=4,
113118
ntapers=2 * nw - 1, taper_weights = fill(1/ntapers, ntapers),
114-
onesided::Bool=T <: Real,
119+
onesided::Bool=T <: Real, ftest::Bool = false,
115120
fft_flags=FFTW.MEASURE) where {T}
116121
if onesided && T <: Complex
117122
throw(ArgumentError("cannot compute one-sided FFT of a complex signal"))
@@ -131,7 +136,7 @@ function MTConfig{T}(n_samples; fs=1, nfft=nextpow(2, n_samples), window=nothing
131136
end
132137

133138
return MTConfig{T}(n_samples, nfft, ntapers, freq, fs, plan, fft_input_tmp,
134-
fft_output_tmp, window, onesided, r)
139+
fft_output_tmp, window, onesided, ftest, r)
135140
end
136141

137142
function allocate_output(config::MTConfig{T}) where {T}
@@ -222,6 +227,27 @@ See also [`mt_pgram`](@ref) and [`MTConfig`](@ref).
222227
"""
223228
mt_pgram!
224229

230+
function ftest(J::Array{ComplexF64,2}, tapers::Matrix)
231+
K = size(tapers,2)
232+
Kodd = 1:2:K
233+
Keven = 2:2:K
234+
235+
Jp = J[:,Kodd]
236+
H0 = sum(tapers[:,Kodd], dims=1)
237+
H0sq = sum(abs2, H0)
238+
JpH0 = dropdims(sum(Jp.*H0, dims=2), dims=2)
239+
A = JpH0 ./ H0sq
240+
Jhat = A * H0
241+
242+
num = (K-1) * abs2.(A) .* H0sq
243+
den = dropdims(sum(abs2.(Jp.-Jhat), dims=2), dims=2) +
244+
dropdims(sum(abs2.(J[:,Keven]), dims=2), dims=2)
245+
Fval = num ./ den
246+
247+
map!(x -> 1.0 - cdf(FDist(2, 2*Int64(K)-2), x), Fval, Fval)
248+
return Fval
249+
end
250+
225251
function mt_pgram!(output::AbstractVector, signal::AbstractVector, config::MTConfig)
226252
if length(output) != length(config.freq)
227253
throw(DimensionMismatch(lazy"""Expected `output` to be of length `length(config.freq)`;
@@ -233,12 +259,20 @@ function mt_pgram!(output::AbstractVector, signal::AbstractVector, config::MTCon
233259
end
234260
fft_output = config.fft_output_tmp
235261

262+
config.ftest && (ftest_scratch = similar(fft_output, (size(fft_output)..., config.ntapers)))
263+
236264
output .= 0
237265
for taper_index in 1:(config.ntapers)
238266
mt_fft_tapered!(fft_output, signal, taper_index, config)
267+
config.ftest && (ftest_scratch[:,taper_index] .= fft_output)
239268
fft2pow!(output, fft_output, config.nfft, config.r[taper_index], config.onesided)
240269
end
241-
return Periodogram(output, config.freq)
270+
if config.ftest
271+
Fval = ftest(ftest_scratch, config.window)
272+
return PeriodogramF(output, config.freq, Fval)
273+
else
274+
return Periodogram(output, config.freq)
275+
end
242276
end
243277

244278
#####

‎src/periodograms.jl‎

Lines changed: 18 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ using ..Util, ..Windows
66
using Statistics: mean!
77
export arraysplit, nextfastfft, periodogram,
88
WelchConfig, welch_pgram, welch_pgram!,
9-
spectrogram, power, freq, stft,
9+
spectrogram, power, freq, Fpval, stft,
1010
MTConfig, mt_pgram, mt_pgram!,
1111
MTSpectrogramConfig, mt_spectrogram, mt_spectrogram!,
1212
MTCrossSpectraConfig, mt_cross_power_spectra, mt_cross_power_spectra!,
@@ -287,6 +287,23 @@ struct Periodogram2{T,F1<:Union{Frequencies,AbstractRange},F2<:Union{Frequencies
287287
freq2::F2
288288
end
289289

290+
"""
291+
PeriodogramF{T, F<:Union{Frequencies,AbstractRange}, F2<:Union{Frequencies,AbstractRange}, M<:AbstractMatrix{T}} <: TFR{T}
292+
293+
A PeriodogramF object with fields:
294+
- `power::M`
295+
- `freq1::F`
296+
- `Ftest::P`
297+
See [`power`](@ref), [`freq`](@ref), and [`Fpval`](@ref) for further details.
298+
"""
299+
struct PeriodogramF{T,F<:Union{Frequencies,AbstractRange},P<:Vector{Float64},V<:AbstractVector{T}} <: TFR{T}
300+
power::V
301+
freq::F
302+
Fpval::P
303+
end
304+
305+
Fpval(p::PeriodogramF) = p.Fpval
306+
290307
"""
291308
power(p)
292309

‎test/multitaper.jl‎

Lines changed: 12 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ const epsilon = 10^-3
2424
ntapers = 5
2525
fs = 1
2626
window = rand(n_samples, ntapers)
27+
ftest = false
2728
for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
2829
fft_input_tmp = Vector{T}(undef, nfft)
2930
onesided = T <: Real
@@ -35,42 +36,42 @@ const epsilon = 10^-3
3536
plan_fft(fft_input_tmp; flags=fft_flags)
3637
# Test that the current configuration is valid so we know if it errors later
3738
# it's because we changed it, not that it was always broken
38-
@test MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r) isa MTConfig
39+
@test MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r) isa MTConfig
3940
# now do a series of changes (in let-blocks to introduce new local bindings)
4041
# and check that they are each invalid
4142
let n_samples = 201
42-
@test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
43+
@test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
4344
end
4445
let n_samples = -1
45-
@test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
46+
@test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
4647
end
4748
let ntapers = -1
48-
@test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
49+
@test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
4950
end
5051
let fs = -1
51-
@test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
52+
@test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
5253
end
5354
let fft_input_tmp = Vector{T}(undef, 2*nfft)
54-
@test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
55+
@test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
5556
end
5657
let nfft = 2*nfft, fft_input_tmp = Vector{T}(undef, 2*nfft)
57-
@test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
58+
@test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
5859
end
5960
let window = rand(2*ntapers, n_samples)
60-
@test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
61+
@test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
6162
end
6263
let r = 1.0
63-
@test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
64+
@test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
6465
end
6566
if T <: Complex
6667
let onesided=true, freqs=rfftfreq(nfft, fs), fft_output_tmp=Vector{fftouttype(T)}(undef, length(freqs)), plan=plan_fft(fft_input_tmp; flags=fft_flags)
67-
@test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
68+
@test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
6869
end
6970
end
7071
# let's check that our original config is still valid, meaning we haven't
7172
# accidentally broken the config by not properly scoping our changes
7273
# (and therefore invalidating the tests that follow).
73-
@test MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r) isa MTConfig
74+
@test MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r) isa MTConfig
7475
end
7576
end
7677

0 commit comments

Comments
 (0)