diff --git a/Project.toml b/Project.toml index 34dfc82..ebabe09 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [compat] FFTW = "1" +LinearAlgebra = "1" julia = "1.6" [extras] diff --git a/README.md b/README.md index 7242989..599477e 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ This is a fork of https://github.com/HolyLab/RFFT.jl with a new UUID such that i ---- -In-place real FFTs for Julia +In-place real FFTs for Julia. Supports "plans" to optimize the algorithm for transformations that you'll perform many times. For example ```julia @@ -14,7 +14,14 @@ a = rand(Float64, 100, 150) # initialize a buffer 'RCpair' that contains a real and complex space buf = RFFT.RCpair{Float64}(undef, size(a)) +``` + +`real(buf)` views the underlying memory buffer as an array reals, while `complex(buf)` views the same +memory buffer as an array of complexes. The user is responsible for keeping track of which view is currently relevant. + +If you'll be performing lots of FFTs on this buffer, it's best to create an optimized plan. +```julia # create the plan plan = RFFT.plan_rfft!(buf; flags=FFTW.MEASURE) @@ -24,3 +31,5 @@ copy!(buf, new) new_fft = plan(buf) ``` + +`RCpair` can be used to implement fast convolutions and many other fourier-based operations on real-valued data. diff --git a/src/RFFT.jl b/src/RFFT.jl index 93eee73..8a0715b 100644 --- a/src/RFFT.jl +++ b/src/RFFT.jl @@ -1,3 +1,14 @@ +""" +# RFFT + +Highlights of the RFFT package: +- [`RCpair`](@ref): a buffer for in-place real-to-complex FFTs +- [`plan_rfft!`](@ref): create a plan for in-place real-to-complex FFTs +- [`plan_irfft!`](@ref): create a plan for in-place complex-to-real FFTs + +Use `real(RC)` and `complex(RC)` to access the real and complex views of the buffer, respectively. +`copy!(RC, A)` can be used to copy data (either real- or complex-valued) into the buffer. +""" module RFFT using FFTW, LinearAlgebra @@ -6,12 +17,34 @@ export RCpair, plan_rfft!, plan_irfft!, rfft!, irfft!, normalization import Base: real, complex, copy, copy! -mutable struct RCpair{T<:AbstractFloat,N,RType<:AbstractArray{T,N},CType<:AbstractArray{Complex{T},N}} +struct RCpair{T<:AbstractFloat,N,RType<:AbstractArray{T,N},CType<:AbstractArray{Complex{T},N}} R::RType C::CType dims::Vector{Int} end +""" + RCPair{T<:AbstractFloat}(undef, realsize::Dims, dims=1:length(realsize)) + +Create a buffer for performing in-place fourier transforms (and inverse) on real-valued data. +A single underlying buffer can be viewed either as real data or as complex data: + +```julia +RC = RCpair{Float64}(undef, (10, 10)) +real(RC) # 10×10 real array +complex(RC) # 6×10 complex array +``` + +`dims` can be used to control which dimensions are transformed. + +The user is responsible to keep track of the current state of the buffer: + +```julia +copy!(RC, rand(10, 10)) # copies real-valued data into the buffer; `real(RC)` is the relevant view +rfft!(RC) # computes the FFT of the real-valued data; now `complex(RC)` is the relevant view +irfft!(RC) # computes the inverse FFT of the complex-valued data; now `real(RC)` is the relevant view +``` +""" function RCpair{T}(::UndefInitializer, realsize::Dims{N}, dims=1:length(realsize)) where {T<:AbstractFloat,N} sz = [realsize...] firstdim = dims[1] @@ -41,6 +74,16 @@ rplan_fwd(R, C, dims, flags, tlim) = FFTW.rFFTWPlan{eltype(R),FFTW.FORWARD,true,ndims(R)}(R, C, dims, flags, tlim) rplan_inv(R, C, dims, flags, tlim) = FFTW.rFFTWPlan{eltype(R),FFTW.BACKWARD,true,ndims(R)}(R, C, dims, flags, tlim) + +""" + plan = plan_rfft!(RC::RCpair; flags=FFTW.ESTIMATE, timelimit=FFTW.NO_TIMELIMIT) + +Create a plan for performing the real-to-complex FFT on the data in `RC`. +Perform the FFT with `plan(RC)`. + +Planning allows you to optimize the performance of (I)FFTs for particular sizes of arrays. +See the FFTW documentation for more information about planning. +""" function plan_rfft!(RC::RCpair{T}; flags::Integer = FFTW.ESTIMATE, timelimit::Real = FFTW.NO_TIMELIMIT) where T p = rplan_fwd(RC.R, RC.C, RC.dims, flags, timelimit) return Z::RCpair -> begin @@ -49,6 +92,16 @@ function plan_rfft!(RC::RCpair{T}; flags::Integer = FFTW.ESTIMATE, timelimit::Re return Z end end + +""" + plan = plan_rfft!(RC::RCpair; flags=FFTW.ESTIMATE, timelimit=FFTW.NO_TIMELIMIT) + +Create a plan for performing the real-to-complex FFT on the data in `RC`. +Perform the FFT with `plan(RC)`. + +Planning allows you to optimize the performance of (I)FFTs for particular sizes of arrays. +See the FFTW documentation for more information about planning. +""" function plan_irfft!(RC::RCpair{T}; flags::Integer = FFTW.ESTIMATE, timelimit::Real = FFTW.NO_TIMELIMIT) where T p = rplan_inv(RC.C, RC.R, RC.dims, flags, timelimit) return Z::RCpair -> begin