diff --git a/Project.toml b/Project.toml
index 354be746..fab9445e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ version = "0.8.3"
 
 [deps]
 Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -22,6 +23,7 @@ OffsetArraysExt = "OffsetArrays"
 [compat]
 Bessels = "0.2"
 DelimitedFiles = "1.6"
+Distributions = "0.25.119"
 FFTW = "1.8"
 IterTools = "1.4"
 LinearAlgebra = "1.6"
diff --git a/src/multitaper.jl b/src/multitaper.jl
index afc32378..082515bb 100644
--- a/src/multitaper.jl
+++ b/src/multitaper.jl
@@ -2,6 +2,8 @@
 ##### Multitapered periodogram
 #####
 
+using Distributions
+
 struct MTConfig{T,R1,F,P,T1,T2,W,R2}
     n_samples::Int
     fs::R1
@@ -13,11 +15,12 @@ struct MTConfig{T,R1,F,P,T1,T2,W,R2}
     fft_output_tmp::T2
     window::W
     onesided::Bool
+    ftest::Bool
     r::R2 # inverse normalization; e.g. equal to `fs*N` for an unwindowed/untapered periodogram of a signal of length `N`
           # 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)
           # can be adjusted to weight the tapers, e.g. by the eigenvalues of the DPSS windows
     function MTConfig{T}(n_samples, nfft, ntapers, freq::F, fs::R1, plan::P, fft_input_tmp::T1,
-                         fft_output_tmp::T2, window::W, onesided, r::R2) where {T,R1,F,P,T1,T2,W,R2}
+                         fft_output_tmp::T2, window::W, onesided, ftest, r::R2) where {T,R1,F,P,T1,T2,W,R2}
         n_samples > 0 || throw(ArgumentError("`n_samples` must be positive"))
         nfft >= n_samples || throw(ArgumentError("Must have `nfft >= n_samples`"))
         ntapers > 0 || throw(ArgumentError("`ntapers` must be positive"))
@@ -44,7 +47,7 @@ struct MTConfig{T,R1,F,P,T1,T2,W,R2}
                                         ntapers, freq, plan,
                                         fft_input_tmp,
                                         fft_output_tmp, window,
-                                        onesided, r)
+                                        onesided, ftest, r)
     end
 end
 
@@ -85,6 +88,7 @@ end
             ntapers = 2 * nw - 1,
             taper_weights = fill(1/ntapers, ntapers),
             onesided::Bool=T<:Real,
+            ftest::Bool=false,
             fft_flags = FFTW.MEASURE)
 
 Creates a config object which holds the configuration state
@@ -107,11 +111,12 @@ as none of the input arguments change.
   The default setting is to simply average them.
 * `onesided`: whether or not to compute a "one-sided" FFT by using that real signal data
   yields conjugate-symmetry in Fourier space.
+* `ftest`:  whether or not to compute the significance of harmonic line compenents as described in Thomson (1982; 10.1109/PROC.1982.12433)
 * `fft_flags`: flags to control how the FFT plan is generated.
 """
 function MTConfig{T}(n_samples; fs=1, nfft=nextpow(2, n_samples), window=nothing, nw=4,
                      ntapers=2 * nw - 1, taper_weights = fill(1/ntapers, ntapers),
-                    onesided::Bool=T <: Real,
+                     onesided::Bool=T <: Real, ftest::Bool = false,
                      fft_flags=FFTW.MEASURE) where {T}
     if onesided && T <: Complex
         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
     end
 
     return MTConfig{T}(n_samples, nfft, ntapers, freq, fs, plan, fft_input_tmp,
-                       fft_output_tmp, window, onesided, r)
+                       fft_output_tmp, window, onesided, ftest, r)
 end
 
 function allocate_output(config::MTConfig{T}) where {T}
@@ -222,6 +227,27 @@ See also [`mt_pgram`](@ref) and [`MTConfig`](@ref).
 """
 mt_pgram!
 
+function ftest(J::Array{ComplexF64,2}, tapers::Matrix)
+  K = size(tapers,2)
+  Kodd = 1:2:K
+  Keven = 2:2:K
+
+  Jp = J[:,Kodd]
+  H0 = sum(tapers[:,Kodd], dims=1)
+  H0sq = sum(abs2, H0)
+  JpH0 = dropdims(sum(Jp.*H0, dims=2), dims=2)
+  A = JpH0 ./ H0sq
+  Jhat = A * H0
+
+  num = (K-1) * abs2.(A) .* H0sq
+  den = dropdims(sum(abs2.(Jp.-Jhat), dims=2), dims=2) +
+        dropdims(sum(abs2.(J[:,Keven]), dims=2), dims=2)
+  Fval = num ./ den
+
+  map!(x -> 1.0 - cdf(FDist(2, 2*Int64(K)-2), x), Fval, Fval)
+  return Fval
+end
+
 function mt_pgram!(output::AbstractVector, signal::AbstractVector, config::MTConfig)
     if length(output) != length(config.freq)
         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
     end
     fft_output = config.fft_output_tmp
 
+    config.ftest && (ftest_scratch = similar(fft_output, (size(fft_output)..., config.ntapers)))
+
     output .= 0
     for taper_index in 1:(config.ntapers)
         mt_fft_tapered!(fft_output, signal, taper_index, config)
+        config.ftest && (ftest_scratch[:,taper_index] .= fft_output)
         fft2pow!(output, fft_output, config.nfft, config.r[taper_index], config.onesided)
     end
-    return Periodogram(output, config.freq)
+    if config.ftest
+        Fval = ftest(ftest_scratch, config.window)
+        return PeriodogramF(output, config.freq, Fval)
+    else
+        return Periodogram(output, config.freq)
+    end
 end
 
 #####
diff --git a/src/periodograms.jl b/src/periodograms.jl
index d4dc2a3d..e1ec705e 100644
--- a/src/periodograms.jl
+++ b/src/periodograms.jl
@@ -6,7 +6,7 @@ using ..Util, ..Windows
 using Statistics: mean!
 export arraysplit, nextfastfft, periodogram,
        WelchConfig, welch_pgram, welch_pgram!,
-       spectrogram, power, freq, stft,
+       spectrogram, power, freq, Fpval, stft,
        MTConfig, mt_pgram, mt_pgram!,
        MTSpectrogramConfig, mt_spectrogram, mt_spectrogram!,
        MTCrossSpectraConfig, mt_cross_power_spectra, mt_cross_power_spectra!,
@@ -287,6 +287,23 @@ struct Periodogram2{T,F1<:Union{Frequencies,AbstractRange},F2<:Union{Frequencies
     freq2::F2
 end
 
+"""
+    PeriodogramF{T, F<:Union{Frequencies,AbstractRange}, F2<:Union{Frequencies,AbstractRange}, M<:AbstractMatrix{T}} <: TFR{T}
+
+A PeriodogramF object with fields:
+- `power::M`
+- `freq1::F`
+- `Ftest::P`
+See [`power`](@ref), [`freq`](@ref), and [`Fpval`](@ref) for further details.
+"""
+struct PeriodogramF{T,F<:Union{Frequencies,AbstractRange},P<:Vector{Float64},V<:AbstractVector{T}} <: TFR{T}
+    power::V
+    freq::F
+    Fpval::P
+end
+
+Fpval(p::PeriodogramF) = p.Fpval
+
 """
     power(p)
 
diff --git a/test/multitaper.jl b/test/multitaper.jl
index c7ccb274..5fc16ab7 100644
--- a/test/multitaper.jl
+++ b/test/multitaper.jl
@@ -24,6 +24,7 @@ const epsilon = 10^-3
         ntapers = 5
         fs = 1
         window = rand(n_samples, ntapers)
+        ftest = false
         for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
             fft_input_tmp = Vector{T}(undef, nfft)
             onesided = T <: Real
@@ -35,42 +36,42 @@ const epsilon = 10^-3
                 plan_fft(fft_input_tmp; flags=fft_flags)
             # Test that the current configuration is valid so we know if it errors later
             # it's because we changed it, not that it was always broken
-            @test MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r) isa MTConfig
+            @test MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r) isa MTConfig
             # now do a series of changes (in let-blocks to introduce new local bindings)
             # and check that they are each invalid
             let n_samples = 201
-                @test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
+                @test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
             end
             let n_samples = -1
-                @test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
+                @test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
             end
             let ntapers = -1
-                @test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
+                @test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
             end
             let fs = -1
-                @test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
+                @test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
             end
             let fft_input_tmp = Vector{T}(undef, 2*nfft)
-                @test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
+                @test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
             end
             let nfft = 2*nfft, fft_input_tmp = Vector{T}(undef, 2*nfft)
-                @test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
+                @test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
             end
             let window = rand(2*ntapers, n_samples)
-                @test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
+                @test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
             end
             let r = 1.0
-                @test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
+                @test_throws DimensionMismatch MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
             end
             if T <: Complex
                 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)
-                    @test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r)
+                    @test_throws ArgumentError MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r)
                 end
             end
             # let's check that our original config is still valid, meaning we haven't
             # accidentally broken the config by not properly scoping our changes
             # (and therefore invalidating the tests that follow).
-            @test MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, r) isa MTConfig
+            @test MTConfig{T}(n_samples, nfft, ntapers, freqs, fs, plan, fft_input_tmp, fft_output_tmp, window, onesided, ftest, r) isa MTConfig
         end
     end