FFT + Measurements

Specific Domains Numerics
1

I thought it would be cool if I could calculate the FFT using uncertain data, so I tried using Measurements together with FFTW. I did not get it working.

Instead, I wrote a simple FT function which can be used with Measurements.

I’m not sure if it is mathematical sound as Measurements only seems to work on Floats and FFTs often produces complex numbers.

Maybe someone with more experience in FFTs can check that.
I think this could be useful for signal/image processing.

Here is some example code:

using FFTW
using Test
using Measurements

function fft2(v::AbstractArray)
    k = similar(v,Complex{eltype(v)})
    f = fftfreq(length(v))
    for i in 1:length(k)
        k[i] = sum(v .* exp.((-1im * 2*π.*(i-1)) .*f))
    end
    k
end

a = range(0,10,length = 20)

reference = fft(a)
sample = fft2(a)
am = measurement.(a,1.0)
sample2 = fft2(am)

@test length(sample) == length(reference)
@test all(sample .≈ reference)
@test all(Measurements.value.(sample2) .≈ reference)
3

FFTW.jl only works with standard numerical types because it’s a wrapper around a C library. FastTransforms.jl has an implementation of FFT for any custom Julia type:

 using Measurements, FastTransforms

julia> V = measurement.(randn(10), rand(10) ./ 10)
10-element Array{Measurement{Float64},1}:
    0.49 ± 0.081
   0.967 ± 0.069
    0.21 ± 0.034
   0.256 ± 0.065
   0.232 ± 0.055
   0.764 ± 0.02
  -1.065 ± 0.041
 -0.7481 ± 0.0099
  -0.653 ± 0.081
   0.835 ± 0.013

julia> ifft(fft(V))
10-element Array{Complex{Measurement{Float64}},1}:
     0.49 ± 0.081 + 0.0 ± 2.6e-17im
    0.967 ± 0.069 - 2.66e-16 ± 1.9e-17im
     0.21 ± 0.034 + 7.55e-16 ± 2.9e-17im
    0.256 ± 0.065 + 1.0e-16 ± 2.8e-17im
    0.232 ± 0.055 - 1.55e-16 ± 2.7e-17im
     0.764 ± 0.02 - 2.23e-16 ± 3.2e-17im
   -1.065 ± 0.041 - 1.33e-16 ± 2.5e-17im
 -0.7481 ± 0.0099 + 9.33e-16 ± 4.0e-17im
   -0.653 ± 0.081 - 8.9e-17 ± 3.0e-17im
    0.835 ± 0.013 + 7.55e-16 ± 3.6e-17im

Performance is not great, though

Measurements work with complex numbers, if this is what you mean:

julia> using Measurements

julia> z = complex(1 ± 0.1, 2 ± 0.1)
(1.0 ± 0.1) + (2.0 ± 0.1)im

julia> z ^ 2
(-3.0 ± 0.45) + (4.0 ± 0.45)im

julia> sqrt(z ^ 2)
(1.0 ± 0.1) + (2.0 ± 0.1)im
5 Likes