diff --git a/Project.toml b/Project.toml
index 84446d5..52c87a9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,7 +10,10 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [compat]
-julia = "1.6"
+julia = "1.8"
+Distributions = "0.25.108"
+KernelFunctions = "0.10.63"
+Reexport = "1.2.2"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/examples/approx-prior-sample.jl b/examples/approx-prior-sample.jl
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/examples/approx-prior-sample.jl
@@ -0,0 +1 @@
+
diff --git a/examples/densities.jl b/examples/densities.jl
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/examples/densities.jl
@@ -0,0 +1 @@
+
diff --git a/examples/feature-functions.jl b/examples/feature-functions.jl
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/examples/feature-functions.jl
@@ -0,0 +1 @@
+
diff --git a/src/KernelSpectralDensities.jl b/src/KernelSpectralDensities.jl
index 6bed9c4..06bbaa4 100644
--- a/src/KernelSpectralDensities.jl
+++ b/src/KernelSpectralDensities.jl
@@ -12,11 +12,9 @@ export SpectralDensity
 export ShiftedRFF, DoubleRFF
 export ApproximateGPSample
 
-# write tests to verify spectral density via fourier transforms
-# also add SpectralKernel (which can then be learned?)
-
 include("base.jl")
 include("expkernels.jl")
+include("matern.jl")
 include("features.jl")
 include("approx_prior.jl")
 
diff --git a/src/base.jl b/src/base.jl
index 2901802..9500072 100644
--- a/src/base.jl
+++ b/src/base.jl
@@ -1,8 +1,6 @@
 
 abstract type AbstractSpectralDensity end
 
-(S::AbstractSpectralDensity)(w) = error("Not implemented")
-
 """
     rand(S::AbstractSpectralDensity, [n::Int])
 
@@ -37,22 +35,63 @@ julia> k = SqExponentialKernel();
 julia> S = SpectralDensity(k, 1);
 
 julia> S(0.0)
-2.5066282746310002
+2.5066282746310007
 
 julia> S = SpectralDensity(k, 2);
 
 julia> S(zeros(2))
-6.283185307179585
+6.2831853071795845
 ```
 """
-struct SpectralDensity{K<:KernelFunctions.Kernel} <: AbstractSpectralDensity
+struct SpectralDensity{K<:KernelFunctions.Kernel,D<:Distribution} <: AbstractSpectralDensity
     kernel::K
-    dim::Int
+    # dim::Int
+    d::D
 
     function SpectralDensity(kernel::KernelFunctions.Kernel, dim::Int)
         if dim < 1
             throw(ArgumentError("Dimension must be greater than 0"))
         end
-        return new{typeof(kernel)}(kernel, dim)
+
+        sk, l = _deconstruct_kernel(kernel, dim)
+        d = _spectral_distribution(sk, l)
+
+        return new{typeof(kernel),typeof(d)}(kernel, d)
+    end
+end
+
+function (S::SpectralDensity)(w)
+    return pdf(S.d, w)
+end
+
+function rand(rng::AbstractRNG, S::SpectralDensity, n::Int...)
+    return rand(rng, S.d, n...)
+end
+
+# ToDo: This could perhaps go into a separate file
+function _deconstruct_kernel(ker::KernelFunctions.SimpleKernel, dim::Int)
+    if dim == 1
+        l = 1.0
+    else
+        l = ones(dim)
+    end
+    return ker, l
+end
+
+function _deconstruct_kernel(
+    ker::TransformedKernel{<:KernelFunctions.SimpleKernel,<:ScaleTransform}, dim::Int
+)
+    l = inv(only(ker.transform.s))
+    if dim > 1
+        l = ones(dim) * l
     end
+    return ker.kernel, l
+end
+
+function _deconstruct_kernel(ker::TransformedKernel, dim::Int)
+    return throw(MethodError(_deconstruct_kernel, (ker, dim)))
 end
+
+function _spectral_distribution(ker::KernelFunctions.Kernel, l)
+    return throw(MethodError(_spectral_distribution, (ker,)))
+end
\ No newline at end of file
diff --git a/src/expkernels.jl b/src/expkernels.jl
index 09b277c..c78aa20 100644
--- a/src/expkernels.jl
+++ b/src/expkernels.jl
@@ -3,40 +3,11 @@
 ## Squared ExponentialKernel
 
 # ToDo: Not sure about distances? Do all work?
-(S::SpectralDensity{<:SqExponentialKernel})(w) = _sqexp(w, 1)
-
-function (S::SpectralDensity{<:TransformedKernel{<:SqExponentialKernel,<:ScaleTransform}})(
-    w
-)
-    l = 1 / only(S.kernel.transform.s)
-    return _sqexp(w, l^2)
-end
-
-_sqexp(w::Real, l2::Real) = sqrt(2 * l2 * π) * exp(-2 * l2 * π^2 * w^2)
-function _sqexp(w::AbstractVector{<:Real}, l2::Real)
-    d = length(w)
-    return sqrt(2 * l2 * π)^d * exp(-2 * l2 * π^2 * dot(w, w))
-end
-
-function rand(rng::AbstractRNG, S::SpectralDensity{<:SqExponentialKernel}, n::Int...)
-    return _sqexprand(rng, S.dim, 1, n...)
-end
-
-function rand(
-    rng::AbstractRNG,
-    S::SpectralDensity{<:TransformedKernel{<:SqExponentialKernel,<:ScaleTransform}},
-    n::Int...,
-)
-    l = 1 / only(S.kernel.transform.s)
-    return _sqexprand(rng, S.dim, l, n...)
+function _spectral_distribution(ker::SqExponentialKernel, l::Real)
+    return inv(2 * π * l) * Normal()
 end
 
-function _sqexprand(rng::AbstractRNG, d::Int, l::Real, n::Int...)
-    σ = 1 / (2 * l * π)
-    if d == 1
-        return rand(rng, Normal(0, σ), n...)
-    elseif d > 1
-        σv = ones(d) * abs2(σ)
-        return rand(rng, MvNormal(Diagonal(σv)), n...)
-    end
+function _spectral_distribution(ker::SqExponentialKernel, l::AbstractVector)
+    σv = abs2.(inv.(2 * π * l))
+    return MvNormal(Diagonal(σv))
 end
diff --git a/src/matern.jl b/src/matern.jl
new file mode 100644
index 0000000..1c0024f
--- /dev/null
+++ b/src/matern.jl
@@ -0,0 +1,24 @@
+
+################################################### 
+## Matern Kernels
+
+MaternKernels = Union{MaternKernel,Matern32Kernel,Matern52Kernel}
+
+_matern_order(k::MaternKernel) = only(k.ν)
+_matern_order(::Matern32Kernel) = 3 / 2
+_matern_order(::Matern52Kernel) = 5 / 2
+
+# rewrite everything as returning a distribution (kind of as originally planned)
+# should be able to abstract/ generalize a lot of the special casing
+function _spectral_distribution(kernel::MaternKernels, l::Real)
+    ν = _matern_order(kernel)
+    return inv(2 * π * l) * TDist(2 * ν)
+end
+
+function _spectral_distribution(kernel::MaternKernels, l::AbstractVector)
+    ν = _matern_order(kernel)
+    n = length(l)
+    l = inv.(2 * π * l) .^ 2
+    D = Distributions.MvTDist(2 * ν, zeros(n), diagm(l))
+    return D
+end
diff --git a/test/base.jl b/test/base.jl
index e120726..a63de88 100644
--- a/test/base.jl
+++ b/test/base.jl
@@ -2,11 +2,13 @@ using KernelSpectralDensities
 using Test
 
 @testset "Fallback" begin
-    ker = ZeroKernel()
+    ker = ConstantKernel()
+
+    @test_throws MethodError SpectralDensity(ker, 1)
 
-    S = SpectralDensity(ker, 1)
+    kert = ker ∘ SelectTransform(1.0)
 
-    @test_throws ErrorException S(1.0)
+    @test_throws MethodError SpectralDensity(kert, 1)
 end
 
 @testset "Dimension check" begin
diff --git a/test/expkernels.jl b/test/expkernels.jl
index c136de1..f1b7428 100644
--- a/test/expkernels.jl
+++ b/test/expkernels.jl
@@ -1,36 +1,48 @@
+if (!@isdefined RUN_TESTS) || !RUN_TESTS
+    using CairoMakie
+    show_plot = true
+    include("test_utils.jl")
+else
+    show_plot = false
+end
+
 @testset "SquaredExponential Kernel" begin
     ker = SqExponentialKernel()
     @testset "1D" begin
         @testset "Pure" begin
-            w_interval = [-2.0, 2.0]
+            # ker = SqExponentialKernel()
+            w_interval = 2.0
             t_interval = [0.0, 4.0]
 
-            f = test_spectral_density(ker, w_interval, t_interval)
+            test_spectral_density(ker, w_interval, t_interval; show_plot)
         end
 
         @testset "Scaled" begin
-            ker = with_lengthscale(SqExponentialKernel(), 0.7)
-            w_interval = [-2.0, 2.0]
+            # ker = SqExponentialKernel()
+            kert = with_lengthscale(ker, 0.7)
+            w_interval = 2.0
             t_interval = [0.0, 4.0]
 
-            f = test_spectral_density(ker, w_interval, t_interval)
+            f = test_spectral_density(kert, w_interval, t_interval; show_plot)
         end
     end
 
     @testset "2D" begin
         @testset "Pure" begin
+            # ker = SqExponentialKernel()
             w_interval = [-2.0, 2.0]
             x_interval = [-2.0, 2.0]
 
-            f = test_2Dspectral_density(ker, w_interval, x_interval)
+            f = test_2Dspectral_density(ker, w_interval, x_interval; show_plot)
         end
 
         @testset "Scaled" begin
-            ker = with_lengthscale(SqExponentialKernel(), 0.7)
+            # ker = SqExponentialKernel()
+            kert = with_lengthscale(ker, 0.7)
             w_interval = [-2.0, 2.0]
             x_interval = [-2.0, 2.0]
 
-            f = test_2Dspectral_density(ker, w_interval, x_interval)
+            f = test_2Dspectral_density(kert, w_interval, x_interval; show_plot)
         end
     end
 end
diff --git a/test/matern.jl b/test/matern.jl
new file mode 100644
index 0000000..591bbc6
--- /dev/null
+++ b/test/matern.jl
@@ -0,0 +1,128 @@
+if (!@isdefined RUN_TESTS) || !RUN_TESTS
+    using CairoMakie
+    show_plot = true
+    include("test_utils.jl")
+else
+    show_plot = false
+end
+
+@testset "Matern32" begin
+    ker = Matern32Kernel()
+    @testset "1D" begin
+        @testset "Pure" begin
+            w_interval = 1.5
+            t_interval = [0.0, 4.0]
+
+            test_spectral_density(ker, w_interval, t_interval; show_plot)
+        end
+
+        @testset "Scaled" begin
+            # ker = Matern32Kernel()
+            kert = with_lengthscale(ker, 0.7)
+            w_interval = 2.0
+            t_interval = [0.0, 3.5]
+
+            test_spectral_density(kert, w_interval, t_interval; show_plot)
+        end
+    end
+
+    @testset "2D" begin
+        @testset "Pure" begin
+            ker = Matern32Kernel()
+            w_interval = [-1.5, 1.5]
+            x_interval = [-3.0, 3.0]
+
+            test_2Dspectral_density(ker, w_interval, x_interval; show_plot)
+        end
+
+        @testset "Scaled" begin
+            # ker = Matern32Kernel()
+            kert = with_lengthscale(ker, 0.6)
+            w_interval = [-2.0, 2.0]
+            x_interval = [-3.5, 3.5]
+
+            f = test_2Dspectral_density(kert, w_interval, x_interval; show_plot)
+        end
+    end
+end
+
+@testset "Matern52" begin
+    ker = Matern52Kernel()
+    @testset "1D" begin
+        @testset "Pure" begin
+            # ker = Matern52Kernel()
+            w_interval = 1.0
+            t_interval = [0.0, 4.0]
+
+            test_spectral_density(ker, w_interval, t_interval; show_plot)
+        end
+
+        @testset "Scaled" begin
+            # ker = Matern52Kernel()
+            kert = with_lengthscale(ker, 0.7)
+            w_interval = 1.5
+            t_interval = [0.0, 3.5]
+
+            test_spectral_density(kert, w_interval, t_interval; show_plot)
+        end
+    end
+
+    @testset "2D" begin
+        @testset "Pure" begin
+            # ker = Matern52Kernel()
+            w_interval = [-1.1, 1.1]
+            x_interval = [-3.0, 3.0]
+
+            test_2Dspectral_density(ker, w_interval, x_interval; show_plot)
+        end
+
+        @testset "Scaled" begin
+            # ker = Matern52Kernel()
+            kert = with_lengthscale(ker, 0.6)
+            w_interval = [-2.0, 2.0]
+            x_interval = [-3.5, 3.5]
+
+            f = test_2Dspectral_density(kert, w_interval, x_interval; show_plot)
+        end
+    end
+end
+
+@testset "Matern72" begin
+    ker = MaternKernel(; ν=7 / 2)
+    @testset "1D" begin
+        @testset "Pure" begin
+            # ker = MaternKernel(; ν=7 / 2)
+            w_interval = 1.0
+            t_interval = [0.0, 6.0]
+
+            test_spectral_density(ker, w_interval, t_interval; show_plot)
+        end
+
+        @testset "Scaled" begin
+            # ker = MaternKernel(; ν=7 / 2)
+            kert = with_lengthscale(ker, 0.7)
+            w_interval = 1.2
+            t_interval = [0.0, 3.5]
+
+            test_spectral_density(kert, w_interval, t_interval; show_plot)
+        end
+    end
+
+    @testset "2D" begin
+        @testset "Pure" begin
+            # ker = MaternKernel(; ν=7 / 2)
+            w_interval = [-1.1, 1.1]
+            x_interval = [-3.0, 3.0]
+
+            test_2Dspectral_density(ker, w_interval, x_interval; show_plot)
+        end
+    end
+    @testset "Scaled" begin
+        # ker = MaternKernel(; ν=7 / 2)
+        kert = with_lengthscale(ker, 0.6)
+        w_interval = [-2.0, 2.0]
+        x_interval = [-3.5, 3.5]
+
+        f = test_2Dspectral_density(kert, w_interval, x_interval; show_plot)
+    end
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index ddc2eee..62edaf1 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,6 +5,8 @@ using FastGaussQuadrature
 using StatsBase
 using StableRNGs
 
+const RUN_TESTS = true
+
 @info "Packages Loaded"
 
 include("test_utils.jl")
@@ -15,6 +17,8 @@ end
 
 @testset "SpectralDensities" begin
     include("expkernels.jl")
+
+    include("matern.jl")
 end
 
 @testset "Feature functions" begin
diff --git a/test/test_utils.jl b/test/test_utils.jl
index 72bfb4c..4d106da 100644
--- a/test/test_utils.jl
+++ b/test/test_utils.jl
@@ -1,4 +1,11 @@
-function test_spectral_density(ker::Kernel, w_interval, t_interval, plot::Bool=false)
+using KernelSpectralDensities
+using Test
+using LinearAlgebra
+using FastGaussQuadrature
+using StatsBase
+using StableRNGs
+
+function test_spectral_density(ker::Kernel, w_max, t_interval; show_plot::Bool=false)
     rng = StableRNG(123)
     S = SpectralDensity(ker, 1)
     k(t) = ker(0, t)
@@ -7,15 +14,16 @@ function test_spectral_density(ker::Kernel, w_interval, t_interval, plot::Bool=f
     τv = range(t_interval...; length=50)
 
     wv, weights = gausslegendre(400)
-    wv = (wv .+ 1) ./ 2 * (w_interval[2] - w_interval[1]) .+ w_interval[1]
-    c = (w_interval[2] - w_interval[1]) / 2
+    wv = (wv .+ 1) ./ 2 * w_max
+    c = w_max / 2
 
-    ks(t) = c * sum(S.(wv) .* cos.(2 * π * wv * t) .* weights)
+    ks(t) = 2 * c * sum(S.(wv) .* cos.(2 * π * wv * t) .* weights)
 
-    @test norm(k.(τv) .- ks.(τv)) < 1e-5
+    @test norm(k.(τv) .- ks.(τv)) < 5e-3
 
     w_samples = rand(rng, S, Int(1e6))
-    h = fit(Histogram, w_samples; nbins=100)
+    w_bins = range(-w_max, w_max; length=100)
+    h = fit(Histogram, w_samples, w_bins)
     h = normalize(h; mode=:pdf)
 
     midpoints = (h.edges[1][1:(end - 1)] .+ h.edges[1][2:end]) ./ 2
@@ -27,7 +35,7 @@ function test_spectral_density(ker::Kernel, w_interval, t_interval, plot::Bool=f
 
     @test norm(k.(τv) .- kss.(τv)) < 0.1
 
-    if plot
+    if show_plot
         f = Figure()
         ax = Axis(f[1, 1])
         lines!(ax, τv, k.(τv); label="kernel")
@@ -35,7 +43,8 @@ function test_spectral_density(ker::Kernel, w_interval, t_interval, plot::Bool=f
         axislegend()
 
         ax2 = Axis(f[1, 2])
-        lines!(ax2, wv, S.(wv); label="spectral density")
+        pwv = vcat(-reverse(wv), wv)
+        lines!(ax2, pwv, S.(pwv); label="spectral density")
         barplot!(
             ax2,
             midpoints,
@@ -45,11 +54,11 @@ function test_spectral_density(ker::Kernel, w_interval, t_interval, plot::Bool=f
             alpha=0.5,
             label="samples",
         )
-        f
+        display(f)
     end
 end
 
-function test_2Dspectral_density(ker::Kernel, w_interval, x_interval; plot::Bool=false)
+function test_2Dspectral_density(ker::Kernel, w_interval, x_interval; show_plot::Bool=false)
     rng = StableRNG(12345)
     S = SpectralDensity(ker, 2)
     k(x) = ker([0, 0], x)
@@ -61,7 +70,7 @@ function test_2Dspectral_density(ker::Kernel, w_interval, x_interval; plot::Bool
 
     ## recover kernel from density
     # Gauss quadrature
-    wv, weights = gausslegendre(200)
+    wv, weights = gausslegendre(250)
     wv = (wv .+ 1) ./ 2 * (w_interval[2] - w_interval[1]) .+ w_interval[1]
     weights = [w_i * w_j for w_i in weights, w_j in weights]
     Wv = [[w_i, w_j] for w_i in wv, w_j in wv]
@@ -70,7 +79,7 @@ function test_2Dspectral_density(ker::Kernel, w_interval, x_interval; plot::Bool
     ks(t) = c^2 * sum(S.(Wv[:]) .* cos.(2 * π * dot.(Wv[:], [t])) .* weights[:])
 
     Ks = [ks(X_i) for X_i in X]
-    @test norm(K .- Ks) < 0.01
+    @test norm(K .- Ks) < 0.015
 
     ## check sampling
     # w = range(w_interval..., length=80)
@@ -82,7 +91,8 @@ function test_2Dspectral_density(ker::Kernel, w_interval, x_interval; plot::Bool
     # w_tmp = (w_tmp[1, :], w_tmp[2, :])
     w_tmp = (w_samples[1, :], w_samples[2, :])
 
-    h = fit(Histogram, w_tmp; nbins=200)
+    w_bins = range(w_interval...; length=100)
+    h = fit(Histogram, w_tmp, (w_bins, w_bins))
     h = normalize(h; mode=:pdf)
     midpoints1 = (h.edges[1][1:(end - 1)] .+ h.edges[1][2:end]) ./ 2
     midpoints2 = (h.edges[2][1:(end - 1)] .+ h.edges[2][2:end]) ./ 2
@@ -99,7 +109,7 @@ function test_2Dspectral_density(ker::Kernel, w_interval, x_interval; plot::Bool
     #ToDo: better would probably be to test if error decreases with more samples
     @test norm(K .- Kss) < 0.3
 
-    if plot
+    if show_plot
         f = Figure(; size=(900, 1000))
         ax = Axis3(f[1, 1])
         contour3d!(
@@ -109,6 +119,7 @@ function test_2Dspectral_density(ker::Kernel, w_interval, x_interval; plot::Bool
             S.(midpoints);
             levels=11,
             color=:orangered3,
+            linewidth=2.5,
             label="spectral density",
         )
         contour3d!(
@@ -118,6 +129,7 @@ function test_2Dspectral_density(ker::Kernel, w_interval, x_interval; plot::Bool
             h.weights;
             levels=11,
             color=:midnightblue,
+            linewidth=2.5,
             linestyle=:dash,
             label="samples",
         )
@@ -132,6 +144,7 @@ function test_2Dspectral_density(ker::Kernel, w_interval, x_interval; plot::Bool
             Ks;
             levels=11,
             color=:midnightblue,
+            linewidth=2.5,
             linestyle=:dash,
             label="spectral approx",
         )
@@ -142,11 +155,12 @@ function test_2Dspectral_density(ker::Kernel, w_interval, x_interval; plot::Bool
             Kss;
             levels=11,
             color=:darkgreen,
+            linewidth=2.5,
             linestyle=:dashdot,
             label="spectral approx (sample)",
         )
         axislegend(ax2)
 
-        f
+        display(f)
     end
 end