From 22a71f20ca4b14b9c20e2960f15e2ca2f344f2d1 Mon Sep 17 00:00:00 2001 From: Alex Arslan Date: Thu, 18 May 2017 16:20:16 -0700 Subject: [PATCH] Remove the FFTW bindings from Base --- .travis.yml | 2 +- DISTRIBUTING.md | 5 +- LICENSE.md | 1 - Make.inc | 17 - Makefile | 5 - README.arm.md | 3 +- README.md | 3 - base/Makefile | 2 - base/deprecated.jl | 43 ++ base/dft.jl | 591 ----------------- base/docs/basedocs.jl | 2 +- base/dsp.jl | 207 ------ base/exports.jl | 37 -- base/fft/FFTW.jl | 802 ----------------------- base/fft/dct.jl | 103 --- base/sysimg.jl | 6 - base/util.jl | 8 - contrib/julia.lang | 1 - contrib/julia.xml | 1 - contrib/mac/macports.make | 1 - contrib/vagrant/Vagrantfile | 4 +- contrib/windows/msys_build.sh | 2 +- deps/Makefile | 11 +- deps/Versions.make | 1 - deps/fftw.mk | 116 ---- doc/src/manual/documentation.md | 2 +- doc/src/manual/noteworthy-differences.md | 2 - doc/src/stdlib/math.md | 51 -- test/choosetests.jl | 4 - test/dsp.jl | 139 ---- test/fft.jl | 352 ---------- test/offsetarray.jl | 25 - 32 files changed, 54 insertions(+), 2495 deletions(-) delete mode 100644 base/dft.jl delete mode 100644 base/dsp.jl delete mode 100644 base/fft/FFTW.jl delete mode 100644 base/fft/dct.jl delete mode 100644 deps/fftw.mk delete mode 100644 test/dsp.jl delete mode 100644 test/fft.jl diff --git a/.travis.yml b/.travis.yml index c6cc7166c9918..2c2638a4ee3fc 100644 --- a/.travis.yml +++ b/.travis.yml @@ -86,7 +86,7 @@ before_install: BUILDOPTS="-j3 USECLANG=1 LLVM_CONFIG=$(brew --prefix llvm39-julia)/bin/llvm-config LLVM_SIZE=$(brew --prefix llvm39-julia)/bin/llvm-size"; BUILDOPTS="$BUILDOPTS VERBOSE=1 USE_BLAS64=0 SUITESPARSE_INC=-I$(brew --prefix suite-sparse-julia)/include FORCE_ASSERTIONS=1"; BUILDOPTS="$BUILDOPTS LIBBLAS=-lopenblas LIBBLASNAME=libopenblas LIBLAPACK=-lopenblas LIBLAPACKNAME=libopenblas"; - for lib in LLVM SUITESPARSE ARPACK BLAS FFTW LAPACK GMP MPFR PCRE LIBUNWIND; do + for lib in LLVM SUITESPARSE ARPACK BLAS LAPACK GMP MPFR PCRE LIBUNWIND; do BUILDOPTS="$BUILDOPTS USE_SYSTEM_$lib=1"; done; export LDFLAGS="-L$(brew --prefix openblas-julia)/lib -L$(brew --prefix suite-sparse-julia)/lib"; diff --git a/DISTRIBUTING.md b/DISTRIBUTING.md index cdceecb728e30..2d0e6916e2037 100644 --- a/DISTRIBUTING.md +++ b/DISTRIBUTING.md @@ -12,9 +12,8 @@ separated most of the notes by OS. Note that while the code for Julia is [MIT-licensed, with a few exceptions](https://github.com/JuliaLang/julia/blob/master/LICENSE.md), the distribution created by the techniques described herein will be -GPL licensed, as various dependent libraries such as `FFTW` and -`SuiteSparse` are GPL licensed. We do hope to have a -non-GPL distribution of Julia in the future. +GPL licensed, as various dependent libraries such as `SuiteSparse` are +GPL licensed. We do hope to have a non-GPL distribution of Julia in the future. Versioning and Git ------------------ diff --git a/LICENSE.md b/LICENSE.md index 1d812d2a7c300..8904e41b74eb0 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -56,7 +56,6 @@ their own licenses: - [OPENLIBM](https://github.com/JuliaLang/openlibm/blob/master/LICENSE.md) [MIT, BSD-2, ISC] - [OPENSPECFUN](https://github.com/JuliaLang/openspecfun) [MIT, public domain] - [FADDEEVA](http://ab-initio.mit.edu/Faddeeva) [MIT] -- [FFTW](http://fftw.org/doc/License-and-Copyright.html) [GPL2+] - [GMP](http://gmplib.org/manual/Copying.html#Copying) [LGPL3+ or GPL2+] - [LIBGIT2](https://github.com/libgit2/libgit2/blob/development/COPYING) [GPL2+ with unlimited linking exception] - [CURL](https://curl.haxx.se/docs/copyright.html) [MIT/X derivative] diff --git a/Make.inc b/Make.inc index 06c3d17e29bff..fd8d1ddb17904 100644 --- a/Make.inc +++ b/Make.inc @@ -34,7 +34,6 @@ USE_SYSTEM_OPENSPECFUN:=0 USE_SYSTEM_DSFMT:=0 USE_SYSTEM_BLAS:=0 USE_SYSTEM_LAPACK:=0 -USE_SYSTEM_FFTW:=0 USE_SYSTEM_GMP:=0 USE_SYSTEM_MPFR:=0 USE_SYSTEM_ARPACK:=0 @@ -53,8 +52,6 @@ USE_LLVM_SHLIB := 1 ## Settings for various Intel tools # Set to 1 to use MKL USE_INTEL_MKL ?= 0 -# Set to 1 to use MKL FFT -USE_INTEL_MKL_FFT ?= 0 # Set to 1 to use Intel LIBM USE_INTEL_LIBM ?= 0 # Set to 1 to enable profiling with Intel VTune Amplifier @@ -832,14 +829,6 @@ LIBLAPACKNAME := liblapack endif endif -ifeq ($(OS), WINNT) -LIBFFTWNAME := libfftw3 -LIBFFTWFNAME := libfftw3f -else -LIBFFTWNAME := libfftw3_threads -LIBFFTWFNAME := libfftw3f_threads -endif - ifeq ($(USE_SYSTEM_LIBM), 1) LIBM := -lm LIBMNAME := libm @@ -1023,12 +1012,6 @@ LIBBLAS := $(MKL_LDFLAGS) LIBLAPACK := $(MKL_LDFLAGS) endif -ifeq ($(USE_INTEL_MKL_FFT), 1) -USE_SYSTEM_FFTW := 1 -LIBFFTWNAME := libmkl_rt -LIBFFTWFNAME := libmkl_rt -endif - ifeq ($(HAVE_SSP),1) JCPPFLAGS += -DHAVE_SSP=1 ifeq ($(USEGCC),1) diff --git a/Makefile b/Makefile index dda792566108a..b3faa68b2164d 100644 --- a/Makefile +++ b/Makefile @@ -245,11 +245,6 @@ JL_PRIVATE_LIBS := ccalltest ifeq ($(USE_GPL_LIBS), 1) JL_PRIVATE_LIBS += suitesparse_wrapper endif -ifeq ($(USE_SYSTEM_FFTW),0) -ifeq ($(USE_GPL_LIBS), 1) -JL_PRIVATE_LIBS += fftw3 fftw3f fftw3_threads fftw3f_threads -endif -endif ifeq ($(USE_SYSTEM_PCRE),0) JL_PRIVATE_LIBS += pcre endif diff --git a/README.arm.md b/README.arm.md index b4320738be8c8..9de38fdeee94e 100644 --- a/README.arm.md +++ b/README.arm.md @@ -27,7 +27,6 @@ adding the following lines in `Make.user`: override USE_SYSTEM_BLAS=1 override USE_SYSTEM_LAPACK=1 override USE_SYSTEM_LIBM=1 -override USE_SYSTEM_FFTW=1 override USE_SYSTEM_GMP=1 override USE_SYSTEM_MPFR=1 override USE_SYSTEM_ARPACK=1 @@ -36,7 +35,7 @@ override USE_SYSTEM_ARPACK=1 The following command will install all the necessary libraries on Ubuntu: ```` -sudo apt-get install libblas3gf liblapack3gf libarpack2 libfftw3-dev libgmp3-dev \ +sudo apt-get install libblas3gf liblapack3gf libarpack2 libgmp3-dev \ libmpfr-dev libblas-dev liblapack-dev cmake gcc-4.8 \ g++-4.8 gfortran libgfortran3 m4 libedit-dev ```` diff --git a/README.md b/README.md index 6cf204e48a403..9502255f7ba10 100644 --- a/README.md +++ b/README.md @@ -277,7 +277,6 @@ Julia uses the following external libraries, which are automatically downloaded - **[AMOS]** — subroutines for computing Bessel and Airy functions. - **[SuiteSparse]** (>= 4.1) — library of linear algebra routines for sparse matrices. - **[ARPACK]** — collection of subroutines designed to solve large, sparse eigenvalue problems. -- **[FFTW]** (>= 3.3) — library for computing fast Fourier transforms very quickly and efficiently. - **[PCRE]** (>= 10.00) — Perl-compatible regular expressions library. - **[GMP]** (>= 5.0) — GNU multiple precision arithmetic library, needed for `BigInt` support. - **[MPFR]** (>= 3.0) — GNU multiple precision floating point library, needed for arbitrary precision floating point (`BigFloat`) support. @@ -309,7 +308,6 @@ Julia uses the following external libraries, which are automatically downloaded [SuiteSparse]: http://faculty.cse.tamu.edu/davis/suitesparse.html [AMOS]: http://netlib.org/amos [ARPACK]: http://forge.scilab.org/index.php/p/arpack-ng -[FFTW]: http://www.fftw.org [PCRE]: http://www.pcre.org [LLVM]: http://www.llvm.org [FemtoLisp]: https://github.com/JeffBezanson/femtolisp @@ -350,7 +348,6 @@ Add the following to the `Make.user` file: USEICC = 1 USEIFC = 1 USE_INTEL_MKL = 1 - USE_INTEL_MKL_FFT = 1 USE_INTEL_LIBM = 1 It is highly recommended to start with a fresh clone of the Julia repository. diff --git a/base/Makefile b/base/Makefile index 04c97fda88e8d..78091ed580eeb 100644 --- a/base/Makefile +++ b/base/Makefile @@ -53,8 +53,6 @@ ifeq ($(USE_GPL_LIBS), 1) else @echo "const USE_GPL_LIBS = false" >> $@ endif - @echo "const libfftw_name = \"$(LIBFFTWNAME)\"" >> $@ - @echo "const libfftwf_name = \"$(LIBFFTWFNAME)\"" >> $@ @echo "const libllvm_version = \"$$($(LLVM_CONFIG_HOST) --version)\"" >> $@ @echo "const VERSION_STRING = \"$(JULIA_VERSION)\"" >> $@ @echo "const TAGGED_RELEASE_BANNER = \"$(TAGGED_RELEASE_BANNER)\"" >> $@ diff --git a/base/deprecated.jl b/base/deprecated.jl index cedc68e0e2a92..6d941213e85dd 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -1342,6 +1342,49 @@ end @deprecate srand(filename::AbstractString, n::Integer=4) srand(read!(filename, Array{UInt32}(Int(n)))) @deprecate MersenneTwister(filename::AbstractString) srand(MersenneTwister(0), read!(filename, Array{UInt32}(Int(4)))) +# PR #21956 +macro _depfftw(f) + quote + function $(esc(f))(args...; kwargs...) + error($f, " has been moved to the package FFTW.jl.\n", + "Run `Pkg.add(\"FFTW\")` to install FFTW on Julia v0.7 and later, ", + "and then run `using FFTW` to load it.") + end + export $(esc(f)) + end +end +function _fftwie(mod::Module) + s = Symbol(mod) + vars = filter(n -> n !== s, names(mod)) + imports = map(v -> Expr(:import, s, v), vars) + return Expr(:toplevel, imports..., Expr(:export, vars...)) +end + +module DFT + for f in [:bfft, :bfft!, :brfft, :dct, :dct!, :fft, :fft!, :fftshift, :idct, :idct!, + :ifft, :ifft!, :ifftshift, :irfft, :plan_bfft, :plan_bfft!, :plan_brfft, + :plan_dct, :plan_dct!, :plan_fft, :plan_fft!, :plan_idct, :plan_idct!, + :plan_ifft, :plan_ifft!, :plan_irfft, :plan_rfft, :rfft] + @eval Base.@_depfftw($f) + end + module FFTW + for f in [:r2r, :r2r!, :plan_r2r, :plan_r2r!] + @eval Base.@_depfftw($f) + end + end +end +@eval _fftwie($DFT) + +const FFTW = DFT.FFTW +export FFTW + +module DSP + for f in [:conv, :conv2, :deconv, :filt, :filt!, :xcorr] + @eval Base.@_depfftw($f) + end +end +@eval _fftwie($DSP) + # END 0.7 deprecations # BEGIN 1.0 deprecations diff --git a/base/dft.jl b/base/dft.jl deleted file mode 100644 index ab4e7e0269409..0000000000000 --- a/base/dft.jl +++ /dev/null @@ -1,591 +0,0 @@ -# This file is a part of Julia. License is MIT: https://julialang.org/license - -module DFT - -# DFT plan where the inputs are an array of eltype T -abstract type Plan{T} end - -import Base: show, summary, size, ndims, length, eltype, - *, A_mul_B!, inv, \, A_ldiv_B! - -eltype(::Type{Plan{T}}) where {T} = T - -# size(p) should return the size of the input array for p -size(p::Plan, d) = size(p)[d] -ndims(p::Plan) = length(size(p)) -length(p::Plan) = prod(size(p))::Int - -############################################################################## -export fft, ifft, bfft, fft!, ifft!, bfft!, - plan_fft, plan_ifft, plan_bfft, plan_fft!, plan_ifft!, plan_bfft!, - rfft, irfft, brfft, plan_rfft, plan_irfft, plan_brfft - -const FFTWFloat = Union{Float32,Float64} -fftwfloat(x) = _fftwfloat(float(x)) -_fftwfloat(::Type{T}) where {T<:FFTWFloat} = T -_fftwfloat(::Type{Float16}) = Float32 -_fftwfloat(::Type{Complex{T}}) where {T} = Complex{_fftwfloat(T)} -_fftwfloat(::Type{T}) where {T} = error("type $T not supported") -_fftwfloat(x::T) where {T} = _fftwfloat(T)(x) - -complexfloat(x::StridedArray{Complex{<:FFTWFloat}}) = x -realfloat(x::StridedArray{<:FFTWFloat}) = x - -# return an Array, rather than similar(x), to avoid an extra copy for FFTW -# (which only works on StridedArray types). -complexfloat(x::AbstractArray{T}) where {T<:Complex} = copy1(typeof(fftwfloat(zero(T))), x) -complexfloat(x::AbstractArray{T}) where {T<:Real} = copy1(typeof(complex(fftwfloat(zero(T)))), x) - -realfloat(x::AbstractArray{T}) where {T<:Real} = copy1(typeof(fftwfloat(zero(T))), x) - -# copy to a 1-based array, using circular permutation -function copy1(::Type{T}, x) where T - y = Array{T}(map(length, indices(x))) - Base.circcopy!(y, x) -end - -to1(x::AbstractArray) = _to1(indices(x), x) -_to1(::Tuple{Base.OneTo,Vararg{Base.OneTo}}, x) = x -_to1(::Tuple, x) = copy1(eltype(x), x) - -# implementations only need to provide plan_X(x, region) -# for X in (:fft, :bfft, ...): -for f in (:fft, :bfft, :ifft, :fft!, :bfft!, :ifft!, :rfft) - pf = Symbol("plan_", f) - @eval begin - $f(x::AbstractArray) = (y = to1(x); $pf(y) * y) - $f(x::AbstractArray, region) = (y = to1(x); $pf(y, region) * y) - $pf(x::AbstractArray; kws...) = (y = to1(x); $pf(y, 1:ndims(y); kws...)) - end -end - -""" - plan_ifft(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf) - -Same as [`plan_fft`](@ref), but produces a plan that performs inverse transforms -[`ifft`](@ref). -""" -plan_ifft - -""" - plan_ifft!(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf) - -Same as [`plan_ifft`](@ref), but operates in-place on `A`. -""" -plan_ifft! - -""" - plan_bfft!(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf) - -Same as [`plan_bfft`](@ref), but operates in-place on `A`. -""" -plan_bfft! - -""" - plan_bfft(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf) - -Same as [`plan_fft`](@ref), but produces a plan that performs an unnormalized -backwards transform [`bfft`](@ref). -""" -plan_bfft - -""" - plan_fft(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf) - -Pre-plan an optimized FFT along given dimensions (`dims`) of arrays matching the shape and -type of `A`. (The first two arguments have the same meaning as for [`fft`](@ref).) -Returns an object `P` which represents the linear operator computed by the FFT, and which -contains all of the information needed to compute `fft(A, dims)` quickly. - -To apply `P` to an array `A`, use `P * A`; in general, the syntax for applying plans is much -like that of matrices. (A plan can only be applied to arrays of the same size as the `A` -for which the plan was created.) You can also apply a plan with a preallocated output array `Â` -by calling `A_mul_B!(Â, plan, A)`. (For `A_mul_B!`, however, the input array `A` must -be a complex floating-point array like the output `Â`.) You can compute the inverse-transform plan by `inv(P)` -and apply the inverse plan with `P \\ Â` (the inverse plan is cached and reused for -subsequent calls to `inv` or `\\`), and apply the inverse plan to a pre-allocated output -array `A` with `A_ldiv_B!(A, P, Â)`. - -The `flags` argument is a bitwise-or of FFTW planner flags, defaulting to `FFTW.ESTIMATE`. -e.g. passing `FFTW.MEASURE` or `FFTW.PATIENT` will instead spend several seconds (or more) -benchmarking different possible FFT algorithms and picking the fastest one; see the FFTW -manual for more information on planner flags. The optional `timelimit` argument specifies a -rough upper bound on the allowed planning time, in seconds. Passing `FFTW.MEASURE` or -`FFTW.PATIENT` may cause the input array `A` to be overwritten with zeros during plan -creation. - -[`plan_fft!`](@ref) is the same as [`plan_fft`](@ref) but creates a -plan that operates in-place on its argument (which must be an array of complex -floating-point numbers). [`plan_ifft`](@ref) and so on are similar but produce -plans that perform the equivalent of the inverse transforms [`ifft`](@ref) and so on. -""" -plan_fft - -""" - plan_fft!(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf) - -Same as [`plan_fft`](@ref), but operates in-place on `A`. -""" -plan_fft! - -""" - rfft(A [, dims]) - -Multidimensional FFT of a real array `A`, exploiting the fact that the transform has -conjugate symmetry in order to save roughly half the computational time and storage costs -compared with [`fft`](@ref). If `A` has size `(n_1, ..., n_d)`, the result has size -`(div(n_1,2)+1, ..., n_d)`. - -The optional `dims` argument specifies an iterable subset of one or more dimensions of `A` -to transform, similar to [`fft`](@ref). Instead of (roughly) halving the first -dimension of `A` in the result, the `dims[1]` dimension is (roughly) halved in the same way. -""" -rfft - -""" - ifft!(A [, dims]) - -Same as [`ifft`](@ref), but operates in-place on `A`. -""" -ifft! - -""" - ifft(A [, dims]) - -Multidimensional inverse FFT. - -A one-dimensional inverse FFT computes - -```math -\\operatorname{IDFT}(A)[k] = \\frac{1}{\\operatorname{length}(A)} -\\sum_{n=1}^{\\operatorname{length}(A)} \\exp\\left(+i\\frac{2\\pi (n-1)(k-1)} -{\\operatorname{length}(A)} \\right) A[n]. -``` - -A multidimensional inverse FFT simply performs this operation along each transformed dimension of `A`. -""" -ifft - -""" - fft!(A [, dims]) - -Same as [`fft`](@ref), but operates in-place on `A`, which must be an array of -complex floating-point numbers. -""" -fft! - -""" - bfft(A [, dims]) - -Similar to [`ifft`](@ref), but computes an unnormalized inverse (backward) -transform, which must be divided by the product of the sizes of the transformed dimensions -in order to obtain the inverse. (This is slightly more efficient than [`ifft`](@ref) -because it omits a scaling step, which in some applications can be combined with other -computational steps elsewhere.) - -```math -\\operatorname{BDFT}(A)[k] = \\operatorname{length}(A) \\operatorname{IDFT}(A)[k] -``` -""" -bfft - -""" - bfft!(A [, dims]) - -Same as [`bfft`](@ref), but operates in-place on `A`. -""" -bfft! - -# promote to a complex floating-point type (out-of-place only), -# so implementations only need Complex{Float} methods -for f in (:fft, :bfft, :ifft) - pf = Symbol("plan_", f) - @eval begin - $f(x::AbstractArray{<:Real}, region=1:ndims(x)) = $f(complexfloat(x), region) - $pf(x::AbstractArray{<:Real}, region; kws...) = $pf(complexfloat(x), region; kws...) - $f(x::AbstractArray{<:Complex{<:Union{Integer,Rational}}}, region=1:ndims(x)) = $f(complexfloat(x), region) - $pf(x::AbstractArray{<:Complex{<:Union{Integer,Rational}}}, region; kws...) = $pf(complexfloat(x), region; kws...) - end -end -rfft(x::AbstractArray{<:Union{Integer,Rational}}, region=1:ndims(x)) = rfft(realfloat(x), region) -plan_rfft(x::AbstractArray, region; kws...) = plan_rfft(realfloat(x), region; kws...) - -# only require implementation to provide *(::Plan{T}, ::Array{T}) -*(p::Plan{T}, x::AbstractArray) where {T} = p * copy1(T, x) - -# Implementations should also implement A_mul_B!(Y, plan, X) so as to support -# pre-allocated output arrays. We don't define * in terms of A_mul_B! -# generically here, however, because of subtleties for in-place and rfft plans. - -############################################################################## -# To support inv, \, and A_ldiv_B!(y, p, x), we require Plan subtypes -# to have a pinv::Plan field, which caches the inverse plan, and which -# should be initially undefined. They should also implement -# plan_inv(p) to construct the inverse of a plan p. - -# hack from @simonster (in #6193) to compute the return type of plan_inv -# without actually calling it or even constructing the empty arrays. -_pinv_type(p::Plan) = typeof([plan_inv(x) for x in typeof(p)[]]) -pinv_type(p::Plan) = eltype(_pinv_type(p)) - -inv(p::Plan) = - isdefined(p, :pinv) ? p.pinv::pinv_type(p) : (p.pinv = plan_inv(p)) -\(p::Plan, x::AbstractArray) = inv(p) * x -A_ldiv_B!(y::AbstractArray, p::Plan, x::AbstractArray) = A_mul_B!(y, inv(p), x) - -############################################################################## -# implementations only need to provide the unnormalized backwards FFT, -# similar to FFTW, and we do the scaling generically to get the ifft: - -mutable struct ScaledPlan{T,P,N} <: Plan{T} - p::P - scale::N # not T, to avoid unnecessary promotion to Complex - pinv::Plan - ScaledPlan{T,P,N}(p, scale) where {T,P,N} = new(p, scale) -end -ScaledPlan{T}(p::P, scale::N) where {T,P,N} = ScaledPlan{T,P,N}(p, scale) -ScaledPlan(p::Plan{T}, scale::Number) where {T} = ScaledPlan{T}(p, scale) -ScaledPlan(p::ScaledPlan, α::Number) = ScaledPlan(p.p, p.scale * α) - -size(p::ScaledPlan) = size(p.p) - -show(io::IO, p::ScaledPlan) = print(io, p.scale, " * ", p.p) -summary(p::ScaledPlan) = string(p.scale, " * ", summary(p.p)) - -*(p::ScaledPlan, x::AbstractArray) = scale!(p.p * x, p.scale) - -*(α::Number, p::Plan) = ScaledPlan(p, α) -*(p::Plan, α::Number) = ScaledPlan(p, α) -*(I::UniformScaling, p::ScaledPlan) = ScaledPlan(p, I.λ) -*(p::ScaledPlan, I::UniformScaling) = ScaledPlan(p, I.λ) -*(I::UniformScaling, p::Plan) = ScaledPlan(p, I.λ) -*(p::Plan, I::UniformScaling) = ScaledPlan(p, I.λ) - -# Normalization for ifft, given unscaled bfft, is 1/prod(dimensions) -normalization(T, sz, region) = one(T) / Int(prod([sz...][[region...]])) -normalization(X, region) = normalization(real(eltype(X)), size(X), region) - -plan_ifft(x::AbstractArray, region; kws...) = - ScaledPlan(plan_bfft(x, region; kws...), normalization(x, region)) -plan_ifft!(x::AbstractArray, region; kws...) = - ScaledPlan(plan_bfft!(x, region; kws...), normalization(x, region)) - -plan_inv(p::ScaledPlan) = ScaledPlan(plan_inv(p.p), inv(p.scale)) - -A_mul_B!(y::AbstractArray, p::ScaledPlan, x::AbstractArray) = - scale!(p.scale, A_mul_B!(y, p.p, x)) - -############################################################################## -# Real-input DFTs are annoying because the output has a different size -# than the input if we want to gain the full factor-of-two(ish) savings -# For backward real-data transforms, we must specify the original length -# of the first dimension, since there is no reliable way to detect this -# from the data (we can't detect whether the dimension was originally even -# or odd). - -for f in (:brfft, :irfft) - pf = Symbol("plan_", f) - @eval begin - $f(x::AbstractArray, d::Integer) = $pf(x, d) * x - $f(x::AbstractArray, d::Integer, region) = $pf(x, d, region) * x - $pf(x::AbstractArray, d::Integer;kws...) = $pf(x, d, 1:ndims(x);kws...) - end -end - -for f in (:brfft, :irfft) - @eval begin - $f(x::AbstractArray{<:Real}, d::Integer, region=1:ndims(x)) = $f(complexfloat(x), d, region) - $f(x::AbstractArray{<:Complex{<:Union{Integer,Rational}}}, d::Integer, region=1:ndims(x)) = $f(complexfloat(x), d, region) - end -end - -""" - irfft(A, d [, dims]) - -Inverse of [`rfft`](@ref): for a complex array `A`, gives the corresponding real -array whose FFT yields `A` in the first half. As for [`rfft`](@ref), `dims` is an -optional subset of dimensions to transform, defaulting to `1:ndims(A)`. - -`d` is the length of the transformed real array along the `dims[1]` dimension, which must -satisfy `div(d,2)+1 == size(A,dims[1])`. (This parameter cannot be inferred from `size(A)` -since both `2*size(A,dims[1])-2` as well as `2*size(A,dims[1])-1` are valid sizes for the -transformed real array.) -""" -irfft - -""" - brfft(A, d [, dims]) - -Similar to [`irfft`](@ref) but computes an unnormalized inverse transform (similar -to [`bfft`](@ref)), which must be divided by the product of the sizes of the -transformed dimensions (of the real output array) in order to obtain the inverse transform. -""" -brfft - -function rfft_output_size(x::AbstractArray, region) - d1 = first(region) - osize = [size(x)...] - osize[d1] = osize[d1]>>1 + 1 - return osize -end - -function brfft_output_size(x::AbstractArray, d::Integer, region) - d1 = first(region) - osize = [size(x)...] - @assert osize[d1] == d>>1 + 1 - osize[d1] = d - return osize -end - -plan_irfft(x::AbstractArray{Complex{T}}, d::Integer, region; kws...) where {T} = - ScaledPlan(plan_brfft(x, d, region; kws...), - normalization(T, brfft_output_size(x, d, region), region)) - -""" - plan_irfft(A, d [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf) - -Pre-plan an optimized inverse real-input FFT, similar to [`plan_rfft`](@ref) -except for [`irfft`](@ref) and [`brfft`](@ref), respectively. The first -three arguments have the same meaning as for [`irfft`](@ref). -""" -plan_irfft - -############################################################################## - -export fftshift, ifftshift - -fftshift(x) = circshift(x, div.([size(x)...],2)) - -""" - fftshift(x) - -Swap the first and second halves of each dimension of `x`. -""" -fftshift(x) - -function fftshift(x,dim) - s = zeros(Int,ndims(x)) - for i in dim - s[i] = div(size(x,i),2) - end - circshift(x, s) -end - -""" - fftshift(x,dim) - -Swap the first and second halves of the given dimension or iterable of dimensions of array `x`. -""" -fftshift(x,dim) - -ifftshift(x) = circshift(x, div.([size(x)...],-2)) - -""" - ifftshift(x, [dim]) - -Undoes the effect of `fftshift`. -""" -ifftshift - -function ifftshift(x,dim) - s = zeros(Int,ndims(x)) - for i in dim - s[i] = -div(size(x,i),2) - end - circshift(x, s) -end - -############################################################################## - -# FFTW module (may move to an external package at some point): -""" - fft(A [, dims]) - -Performs a multidimensional FFT of the array `A`. The optional `dims` argument specifies an -iterable subset of dimensions (e.g. an integer, range, tuple, or array) to transform along. -Most efficient if the size of `A` along the transformed dimensions is a product of small -primes; see `nextprod()`. See also `plan_fft()` for even greater efficiency. - -A one-dimensional FFT computes the one-dimensional discrete Fourier transform (DFT) as -defined by - -```math -\\operatorname{DFT}(A)[k] = - \\sum_{n=1}^{\\operatorname{length}(A)} - \\exp\\left(-i\\frac{2\\pi - (n-1)(k-1)}{\\operatorname{length}(A)} \\right) A[n]. -``` - -A multidimensional FFT simply performs this operation along each transformed dimension of `A`. - -!!! note - * Julia starts FFTW up with 1 thread by default. Higher performance is usually possible by - increasing number of threads. Use `FFTW.set_num_threads(Sys.CPU_CORES)` to use as many - threads as cores on your system. - - * This performs a multidimensional FFT by default. FFT libraries in other languages such as - Python and Octave perform a one-dimensional FFT along the first non-singleton dimension - of the array. This is worth noting while performing comparisons. For more details, - refer to the [Noteworthy Differences from other Languages](@ref) - section of the manual. -""" -fft - -""" - plan_rfft(A [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf) - -Pre-plan an optimized real-input FFT, similar to [`plan_fft`](@ref) except for -[`rfft`](@ref) instead of [`fft`](@ref). The first two arguments, and the -size of the transformed result, are the same as for [`rfft`](@ref). -""" -plan_rfft - -""" - plan_brfft(A, d [, dims]; flags=FFTW.ESTIMATE; timelimit=Inf) - -Pre-plan an optimized real-input unnormalized transform, similar to -[`plan_rfft`](@ref) except for [`brfft`](@ref) instead of -[`rfft`](@ref). The first two arguments and the size of the transformed result, are -the same as for [`brfft`](@ref). -""" -plan_brfft - -module FFTW - import ..DFT: fft, bfft, ifft, rfft, brfft, irfft, plan_fft, plan_bfft, plan_ifft, - plan_rfft, plan_brfft, plan_irfft, fft!, bfft!, ifft!, plan_fft!, plan_bfft!, plan_ifft!, - Plan, rfft_output_size, brfft_output_size, plan_inv, normalization, ScaledPlan - - export r2r, r2r!, plan_r2r, plan_r2r! - - """ - plan_dct!(A [, dims [, flags [, timelimit]]]) - - Same as [`plan_dct`](@ref), but operates in-place on `A`. - """ - function plan_dct! end - - """ - plan_idct(A [, dims [, flags [, timelimit]]]) - - Pre-plan an optimized inverse discrete cosine transform (DCT), similar to - [`plan_fft`](@ref) except producing a function that computes - [`idct`](@ref). The first two arguments have the same meaning as for - [`idct`](@ref). - """ - function plan_idct end - - """ - plan_dct(A [, dims [, flags [, timelimit]]]) - - Pre-plan an optimized discrete cosine transform (DCT), similar to - [`plan_fft`](@ref) except producing a function that computes - [`dct`](@ref). The first two arguments have the same meaning as for - [`dct`](@ref). - """ - function plan_dct end - - """ - plan_idct!(A [, dims [, flags [, timelimit]]]) - - Same as [`plan_idct`](@ref), but operates in-place on `A`. - """ - function plan_idct! end - - """ - dct(A [, dims]) - - Performs a multidimensional type-II discrete cosine transform (DCT) of the array `A`, using - the unitary normalization of the DCT. The optional `dims` argument specifies an iterable - subset of dimensions (e.g. an integer, range, tuple, or array) to transform along. Most - efficient if the size of `A` along the transformed dimensions is a product of small primes; - see [`nextprod`](@ref). See also [`plan_dct`](@ref) for even greater - efficiency. - """ - function dct end - - """ - idct(A [, dims]) - - Computes the multidimensional inverse discrete cosine transform (DCT) of the array `A` - (technically, a type-III DCT with the unitary normalization). The optional `dims` argument - specifies an iterable subset of dimensions (e.g. an integer, range, tuple, or array) to - transform along. Most efficient if the size of `A` along the transformed dimensions is a - product of small primes; see [`nextprod`](@ref). See also - [`plan_idct`](@ref) for even greater efficiency. - """ - function idct end - - """ - dct!(A [, dims]) - - Same as [`dct!`](@ref), except that it operates in-place on `A`, which must be an - array of real or complex floating-point values. - """ - function dct! end - - """ - idct!(A [, dims]) - - Same as [`idct!`](@ref), but operates in-place on `A`. - """ - function idct! end - - """ - r2r(A, kind [, dims]) - - Performs a multidimensional real-input/real-output (r2r) transform - of type `kind` of the array `A`, as defined in the FFTW manual. - `kind` specifies either a discrete cosine transform of various types - (`FFTW.REDFT00`, `FFTW.REDFT01`, `FFTW.REDFT10`, or - `FFTW.REDFT11`), a discrete sine transform of various types - (`FFTW.RODFT00`, `FFTW.RODFT01`, `FFTW.RODFT10`, or - `FFTW.RODFT11`), a real-input DFT with halfcomplex-format output - (`FFTW.R2HC` and its inverse `FFTW.HC2R`), or a discrete - Hartley transform (`FFTW.DHT`). The `kind` argument may be - an array or tuple in order to specify different transform types - along the different dimensions of `A`; `kind[end]` is used - for any unspecified dimensions. See the FFTW manual for precise - definitions of these transform types, at http://www.fftw.org/doc. - - The optional `dims` argument specifies an iterable subset of - dimensions (e.g. an integer, range, tuple, or array) to transform - along. `kind[i]` is then the transform type for `dims[i]`, - with `kind[end]` being used for `i > length(kind)`. - - See also [`plan_r2r`](@ref) to pre-plan optimized r2r transforms. - """ - function r2r end - - """ - r2r!(A, kind [, dims]) - - Same as [`r2r`](@ref), but operates in-place on `A`, which must be - an array of real or complex floating-point numbers. - """ - function r2r! end - - """ - plan_r2r!(A, kind [, dims [, flags [, timelimit]]]) - - Similar to [`plan_fft`](@ref), but corresponds to [`r2r!`](@ref). - """ - function plan_r2r! end - - """ - plan_r2r(A, kind [, dims [, flags [, timelimit]]]) - - Pre-plan an optimized r2r transform, similar to [`plan_fft`](@ref) - except that the transforms (and the first three arguments) - correspond to [`r2r`](@ref) and [`r2r!`](@ref), respectively. - """ - function plan_r2r end - - (Base.USE_GPL_LIBS || Base.fftw_vendor() == :mkl) && include(joinpath("fft", "FFTW.jl")) -end - -importall .FFTW -export FFTW, dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct! - -############################################################################## - -end diff --git a/base/docs/basedocs.jl b/base/docs/basedocs.jl index 9134bbc3925f7..fb5d06f5ae10c 100644 --- a/base/docs/basedocs.jl +++ b/base/docs/basedocs.jl @@ -17,7 +17,7 @@ as well many great tutorials and learning resources: https://julialang.org/learning/ For help on a specific function or macro, type `?` followed -by its name, e.g. `?fft`, or `?@time`, and press enter. +by its name, e.g. `?cos`, or `?@time`, and press enter. """ kw"help", kw"?", kw"julia" diff --git a/base/dsp.jl b/base/dsp.jl deleted file mode 100644 index 470852e84b561..0000000000000 --- a/base/dsp.jl +++ /dev/null @@ -1,207 +0,0 @@ -# This file is a part of Julia. License is MIT: https://julialang.org/license - -module DSP - -import Base.trailingsize - -export filt, filt!, deconv, conv, conv2, xcorr - -_zerosi(b,a,T) = zeros(promote_type(eltype(b), eltype(a), T), max(length(a), length(b))-1) - -""" - filt(b, a, x, [si]) - -Apply filter described by vectors `a` and `b` to vector `x`, with an optional initial filter -state vector `si` (defaults to zeros). -""" -function filt(b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number}, - x::AbstractArray{T}, si::AbstractArray{S} = _zerosi(b,a,T)) where {T,S} - filt!(Array{promote_type(eltype(b), eltype(a), T, S)}(size(x)), b, a, x, si) -end - -# in-place filtering: returns results in the out argument, which may shadow x -# (and does so by default) - -""" - filt!(out, b, a, x, [si]) - -Same as [`filt`](@ref) but writes the result into the `out` argument, which may -alias the input `x` to modify it in-place. -""" -function filt!(out::AbstractArray, b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number}, - x::AbstractArray{T}, si::AbstractArray{S,N} = _zerosi(b,a,T)) where {T,S,N} - isempty(b) && throw(ArgumentError("filter vector b must be non-empty")) - isempty(a) && throw(ArgumentError("filter vector a must be non-empty")) - a[1] == 0 && throw(ArgumentError("filter vector a[1] must be nonzero")) - if size(x) != size(out) - throw(ArgumentError("output size $(size(out)) must match input size $(size(x))")) - end - - as = length(a) - bs = length(b) - sz = max(as, bs) - silen = sz - 1 - ncols = trailingsize(x,2) - - if size(si, 1) != silen - throw(ArgumentError("initial state vector si must have max(length(a),length(b))-1 rows")) - end - if N > 1 && trailingsize(si,2) != ncols - throw(ArgumentError("initial state vector si must be a vector or have the same number of columns as x")) - end - - size(x,1) == 0 && return out - sz == 1 && return scale!(out, x, b[1]/a[1]) # Simple scaling without memory - - # Filter coefficient normalization - if a[1] != 1 - norml = a[1] - a ./= norml - b ./= norml - end - # Pad the coefficients with zeros if needed - bs 1 ? col : 1] - if as > 1 - _filt_iir!(out, b, a, x, si, col) - else - _filt_fir!(out, b, x, si, col) - end - end - return out -end - -function _filt_iir!(out, b, a, x, si, col) - silen = length(si) - @inbounds for i=1:size(x, 1) - xi = x[i,col] - val = si[1] + b[1]*xi - for j=1:(silen-1) - si[j] = si[j+1] + b[j+1]*xi - a[j+1]*val - end - si[silen] = b[silen+1]*xi - a[silen+1]*val - out[i,col] = val - end -end - -function _filt_fir!(out, b, x, si, col) - silen = length(si) - @inbounds for i=1:size(x, 1) - xi = x[i,col] - val = si[1] + b[1]*xi - for j=1:(silen-1) - si[j] = si[j+1] + b[j+1]*xi - end - si[silen] = b[silen+1]*xi - out[i,col] = val - end -end - -""" - deconv(b,a) -> c - -Construct vector `c` such that `b = conv(a,c) + r`. -Equivalent to polynomial division. -""" -function deconv(b::StridedVector{T}, a::StridedVector{T}) where T - lb = size(b,1) - la = size(a,1) - if lb < la - return [zero(T)] - end - lx = lb-la+1 - x = zeros(T, lx) - x[1] = 1 - filt(b, a, x) -end - -""" - conv(u,v) - -Convolution of two vectors. Uses FFT algorithm. -""" -function conv(u::StridedVector{T}, v::StridedVector{T}) where T<:Base.LinAlg.BlasFloat - nu = length(u) - nv = length(v) - n = nu + nv - 1 - np2 = n > 1024 ? nextprod([2,3,5], n) : nextpow2(n) - upad = [u; zeros(T, np2 - nu)] - vpad = [v; zeros(T, np2 - nv)] - if T <: Real - p = plan_rfft(upad) - y = irfft((p*upad).*(p*vpad), np2) - else - p = plan_fft!(upad) - y = ifft!((p*upad).*(p*vpad)) - end - return y[1:n] -end -conv(u::StridedVector{T}, v::StridedVector{T}) where {T<:Integer} = round.(Int, conv(float(u), float(v))) -conv(u::StridedVector{<:Integer}, v::StridedVector{<:Base.LinAlg.BlasFloat}) = conv(float(u), v) -conv(u::StridedVector{<:Base.LinAlg.BlasFloat}, v::StridedVector{<:Integer}) = conv(u, float(v)) - -""" - conv2(u,v,A) - -2-D convolution of the matrix `A` with the 2-D separable kernel generated by -the vectors `u` and `v`. -Uses 2-D FFT algorithm. -""" -function conv2(u::StridedVector{T}, v::StridedVector{T}, A::StridedMatrix{T}) where T - m = length(u)+size(A,1)-1 - n = length(v)+size(A,2)-1 - B = zeros(T, m, n) - B[1:size(A,1),1:size(A,2)] = A - u = fft([u;zeros(T,m-length(u))]) - v = fft([v;zeros(T,n-length(v))]) - C = ifft(fft(B) .* (u * v.')) - if T <: Real - return real(C) - end - return C -end - -""" - conv2(B,A) - -2-D convolution of the matrix `B` with the matrix `A`. Uses 2-D FFT algorithm. -""" -function conv2(A::StridedMatrix{T}, B::StridedMatrix{T}) where T - sa, sb = size(A), size(B) - At = zeros(T, sa[1]+sb[1]-1, sa[2]+sb[2]-1) - Bt = zeros(T, sa[1]+sb[1]-1, sa[2]+sb[2]-1) - At[1:sa[1], 1:sa[2]] = A - Bt[1:sb[1], 1:sb[2]] = B - p = plan_fft(At) - C = ifft((p*At).*(p*Bt)) - if T <: Real - return real(C) - end - return C -end -conv2(A::StridedMatrix{T}, B::StridedMatrix{T}) where {T<:Integer} = - round.(Int, conv2(float(A), float(B))) -conv2(u::StridedVector{T}, v::StridedVector{T}, A::StridedMatrix{T}) where {T<:Integer} = - round.(Int, conv2(float(u), float(v), float(A))) - -""" - xcorr(u,v) - -Compute the cross-correlation of two vectors. -""" -function xcorr(u, v) - su = size(u,1); sv = size(v,1) - if su < sv - u = [u;zeros(eltype(u),sv-su)] - elseif sv < su - v = [v;zeros(eltype(v),su-sv)] - end - flipdim(conv(flipdim(u, 1), v), 1) -end - -end # module diff --git a/base/exports.jl b/base/exports.jl index 376a9704b6d55..c4395287a279e 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -2,7 +2,6 @@ export # Modules - FFTW, Meta, Operators, Pkg, @@ -862,42 +861,6 @@ export var, varm, -# signal processing - bfft!, - bfft, - brfft, - conv, - conv2, - dct!, - dct, - deconv, - fft!, - fft, - fftshift, - filt, - filt!, - idct!, - idct, - ifft!, - ifft, - ifftshift, - irfft, - plan_bfft!, - plan_bfft, - plan_brfft, - plan_dct!, - plan_dct, - plan_fft!, - plan_fft, - plan_idct!, - plan_idct, - plan_ifft!, - plan_ifft, - plan_irfft, - plan_rfft, - rfft, - xcorr, - # iteration done, next, diff --git a/base/fft/FFTW.jl b/base/fft/FFTW.jl deleted file mode 100644 index 4ef2f701e0c78..0000000000000 --- a/base/fft/FFTW.jl +++ /dev/null @@ -1,802 +0,0 @@ -# This file is a part of Julia. License is MIT: https://julialang.org/license - -import Base: show, *, convert, unsafe_convert, size, strides, ndims, pointer, A_mul_B! - -export export_wisdom, import_wisdom, import_system_wisdom, forget_wisdom, - MEASURE, DESTROY_INPUT, UNALIGNED, CONSERVE_MEMORY, EXHAUSTIVE, - PRESERVE_INPUT, PATIENT, ESTIMATE, WISDOM_ONLY, NO_TIMELIMIT, - R2HC, HC2R, DHT, REDFT00, REDFT01, REDFT10, REDFT11, - RODFT00, RODFT01, RODFT10, RODFT11, - fftwNumber, fftwReal, fftwComplex, flops - -## FFT: Implement fft by calling fftw. - -const libfftw = Base.libfftw_name -const libfftwf = Base.libfftwf_name - -const version = convert(VersionNumber, split(unsafe_string(cglobal( - (:fftw_version,Base.DFT.FFTW.libfftw), UInt8)), ['-', ' '])[2]) - -## Direction of FFT - -const FORWARD = -1 -const BACKWARD = 1 - -## FFTW Flags from fftw3.h - -const MEASURE = UInt32(0) -const DESTROY_INPUT = UInt32(1 << 0) -const UNALIGNED = UInt32(1 << 1) -const CONSERVE_MEMORY = UInt32(1 << 2) -const EXHAUSTIVE = UInt32(1 << 3) # NO_EXHAUSTIVE is default -const PRESERVE_INPUT = UInt32(1 << 4) # cancels DESTROY_INPUT -const PATIENT = UInt32(1 << 5) # IMPATIENT is default -const ESTIMATE = UInt32(1 << 6) -const WISDOM_ONLY = UInt32(1 << 21) -const NO_SIMD = UInt32(1 << 17) # disable SIMD, useful for benchmarking - -## R2R transform kinds - -const R2HC = 0 -const HC2R = 1 -const DHT = 2 -const REDFT00 = 3 -const REDFT01 = 4 -const REDFT10 = 5 -const REDFT11 = 6 -const RODFT00 = 7 -const RODFT01 = 8 -const RODFT10 = 9 -const RODFT11 = 10 - -let k2s = Dict(R2HC => "R2HC", HC2R => "HC2R", DHT => "DHT", REDFT00 => "REDFT00", REDFT01 => "REDFT01", REDFT10 => "REDFT10", REDFT11 => "REDFT11", RODFT00 => "RODFT00", RODFT01 => "RODFT01", RODFT10 => "RODFT10", RODFT11 => "RODFT11") - global kind2string - kind2string(k::Integer) = k2s[Int(k)] -end - -# FFTW floating-point types: - -const fftwNumber = Union{Float64,Float32,Complex128,Complex64} -const fftwReal = Union{Float64,Float32} -const fftwComplex = Union{Complex128,Complex64} -const fftwDouble = Union{Float64,Complex128} -const fftwSingle = Union{Float32,Complex64} -const fftwTypeDouble = Union{Type{Float64},Type{Complex128}} -const fftwTypeSingle = Union{Type{Float32},Type{Complex64}} - -# For ESTIMATE plans, FFTW allows one to pass NULL for the array pointer, -# since it is not written to. Hence, it is convenient to create an -# array-like type that carries a size and a stride like a "real" array -# but which is converted to C_NULL as a pointer. -struct FakeArray{T,N} <: DenseArray{T,N} - sz::NTuple{N,Int} - st::NTuple{N,Int} -end -size(a::FakeArray) = a.sz -strides(a::FakeArray) = a.st -unsafe_convert(::Type{Ptr{T}}, a::FakeArray{T}) where {T} = convert(Ptr{T}, C_NULL) -pointer(a::FakeArray{T}) where {T} = convert(Ptr{T}, C_NULL) -FakeArray{T,N}(::Type{T}, sz::NTuple{N,Int}) = FakeArray{T,N}(sz, colmajorstrides(sz)) -FakeArray{T}(::Type{T}, sz::Int...) = FakeArray(T, sz) -fakesimilar(flags, X, T) = flags & ESTIMATE != 0 ? FakeArray(T, size(X)) : Array{T}(size(X)) -alignment_of(A::FakeArray) = Int32(0) - -## Julia wrappers around FFTW functions - -# _init_() must be called before any FFTW planning routine. -# -- Once FFTW is split into its own module, this can be called -# in the module __init__(), but for now we must call it lazily -# in every routine that might initialize the FFTW planner. -# -- This initializes FFTW's threads support (defaulting to 1 thread). -# If this isn't called before the FFTW planner is created, then -# FFTW's threads algorithms won't be registered or used at all. -# (Previously, we called fftw_cleanup, but this invalidated existing -# plans, causing issue #19892.) -const threads_initialized = Ref(false) -function _init_() - if !threads_initialized[] - stat = ccall((:fftw_init_threads,libfftw), Int32, ()) - statf = ccall((:fftwf_init_threads,libfftwf), Int32, ()) - if stat == 0 || statf == 0 - error("could not initialize FFTW threads") - end - threads_initialized[] = true - end -end - -# Wisdom - -# Import and export wisdom to/from a single file for all precisions, -# which is more user-friendly than requiring the user to call a -# separate routine depending on the fp precision of the plans. This -# requires a bit of trickness since we have to (a) use the libc file -# I/O routines with fftw_export_wisdom_to_file/import_wisdom_from_file -# (b) we need 256 bytes of space padding between the wisdoms to work -# around FFTW's internal file i/o buffering [see the BUFSZ constant in -# FFTW's api/import-wisdom-from-file.c file]. - -function export_wisdom(fname::AbstractString) - _init_() - f = ccall(:fopen, Ptr{Void}, (Cstring,Cstring), fname, :w) - systemerror("could not open wisdom file $fname for writing", f == C_NULL) - ccall((:fftw_export_wisdom_to_file,libfftw), Void, (Ptr{Void},), f) - ccall(:fputs, Int32, (Ptr{UInt8},Ptr{Void}), " "^256, f) # no NUL, hence no Cstring - ccall((:fftwf_export_wisdom_to_file,libfftwf), Void, (Ptr{Void},), f) - ccall(:fclose, Void, (Ptr{Void},), f) -end - -function import_wisdom(fname::AbstractString) - _init_() - f = ccall(:fopen, Ptr{Void}, (Cstring,Cstring), fname, :r) - systemerror("could not open wisdom file $fname for reading", f == C_NULL) - if ccall((:fftw_import_wisdom_from_file,libfftw),Int32,(Ptr{Void},),f)==0|| - ccall((:fftwf_import_wisdom_from_file,libfftwf),Int32,(Ptr{Void},),f)==0 - error("failed to import wisdom from $fname") - end - ccall(:fclose, Void, (Ptr{Void},), f) -end - -function import_system_wisdom() - _init_() - if ccall((:fftw_import_system_wisdom,libfftw), Int32, ()) == 0 || - ccall((:fftwf_import_system_wisdom,libfftwf), Int32, ()) == 0 - error("failed to import system wisdom") - end -end - -function forget_wisdom() - _init_() - ccall((:fftw_forget_wisdom,libfftw), Void, ()) - ccall((:fftwf_forget_wisdom,libfftwf), Void, ()) -end - -# Threads - -function set_num_threads(nthreads::Integer) - _init_() - ccall((:fftw_plan_with_nthreads,libfftw), Void, (Int32,), nthreads) - ccall((:fftwf_plan_with_nthreads,libfftwf), Void, (Int32,), nthreads) -end - -# pointer type for fftw_plan (opaque pointer) - -struct fftw_plan_struct end -const PlanPtr = Ptr{fftw_plan_struct} - -# Planner timelimits - -const NO_TIMELIMIT = -1.0 # from fftw3.h - -function set_timelimit(precision::fftwTypeDouble,seconds) - _init_() - ccall((:fftw_set_timelimit,libfftw), Void, (Float64,), seconds) -end - -function set_timelimit(precision::fftwTypeSingle,seconds) - _init_() - ccall((:fftwf_set_timelimit,libfftwf), Void, (Float64,), seconds) -end - -# Array alignment mod 16: -# FFTW plans may depend on the alignment of the array mod 16 bytes, -# i.e. the address mod 16 of the first element of the array, in order -# to exploit SIMD operations. Julia arrays are, by default, aligned -# to 16-byte boundaries (address mod 16 == 0), but this may not be -# true for data imported from external C code, or for SubArrays. -# Use the undocumented routine fftw_alignment_of to determine the -# alignment of a given pointer modulo whatever FFTW needs; this -# function will be documented in FFTW 3.3.4. - - -if Base.libfftw_name == "libmkl_rt" - alignment_of(A::StridedArray{<:fftwDouble}) = - convert(Int32, convert(Int64, pointer(A)) % 16) - alignment_of(A::StridedArray{<:fftwSingle}) = - convert(Int32, convert(Int64, pointer(A)) % 16) -else - alignment_of{T<:fftwDouble}(A::StridedArray{T}) = - ccall((:fftw_alignment_of, libfftw), Int32, (Ptr{T},), A) - alignment_of{T<:fftwSingle}(A::StridedArray{T}) = - ccall((:fftwf_alignment_of, libfftwf), Int32, (Ptr{T},), A) -end - -# FFTWPlan (low-level) - -# low-level storage of the FFTW plan, along with the information -# needed to determine whether it is applicable. We need to put -# this into a type to support a finalizer on the fftw_plan. -# K is FORWARD/BACKWARD for forward/backward or r2c/c2r plans, respectively. -# For r2r plans, K is a tuple of the transform kinds along each dimension. -abstract type FFTWPlan{T<:fftwNumber,K,inplace} <: Plan{T} end -for P in (:cFFTWPlan, :rFFTWPlan, :r2rFFTWPlan) # complex, r2c/c2r, and r2r - @eval begin - mutable struct $P{T<:fftwNumber,K,inplace,N} <: FFTWPlan{T,K,inplace} - plan::PlanPtr - sz::NTuple{N,Int} # size of array on which plan operates (Int tuple) - osz::NTuple{N,Int} # size of output array (Int tuple) - istride::NTuple{N,Int} # strides of input - ostride::NTuple{N,Int} # strides of output - ialign::Int32 # alignment mod 16 of input - oalign::Int32 # alignment mod 16 of input - flags::UInt32 # planner flags - region::Any # region (iterable) of dims that are transormed - pinv::ScaledPlan - function $P{T,K,inplace,N}(plan::PlanPtr, flags::Integer, R::Any, - X::StridedArray{T,N}, Y::StridedArray) where {T<:fftwNumber,K,inplace,N} - p = new(plan, size(X), size(Y), strides(X), strides(Y), - alignment_of(X), alignment_of(Y), flags, R) - finalizer(p, destroy_plan) - p - end - end - end -end - -size(p::FFTWPlan) = p.sz - -unsafe_convert(::Type{PlanPtr}, p::FFTWPlan) = p.plan - -destroy_plan(plan::FFTWPlan{<:fftwDouble}) = - ccall((:fftw_destroy_plan,libfftw), Void, (PlanPtr,), plan) - -destroy_plan(plan::FFTWPlan{<:fftwSingle}) = - ccall((:fftwf_destroy_plan,libfftwf), Void, (PlanPtr,), plan) - -cost(plan::FFTWPlan{<:fftwDouble}) = - ccall((:fftw_cost,libfftw), Float64, (PlanPtr,), plan) -cost(plan::FFTWPlan{<:fftwSingle}) = - ccall((:fftwf_cost,libfftwf), Float64, (PlanPtr,), plan) - -function arithmetic_ops(plan::FFTWPlan{<:fftwDouble}) - # Change to individual Ref after we can allocate them on stack - ref = Ref{NTuple{3,Float64}}() - ptr = Ptr{Float64}(Base.unsafe_convert(Ptr{NTuple{3,Float64}}, ref)) - ccall((:fftw_flops,libfftw), Void, - (PlanPtr,Ptr{Float64},Ptr{Float64},Ptr{Float64}), - plan, ptr, ptr + 8, ptr + 16) - (round(Int64, ref[][1]), round(Int64, ref[][2]), round(Int64, ref[][3])) -end -function arithmetic_ops(plan::FFTWPlan{<:fftwSingle}) - # Change to individual Ref after we can allocate them on stack - ref = Ref{NTuple{3,Float64}}() - ptr = Ptr{Float64}(Base.unsafe_convert(Ptr{NTuple{3,Float64}}, ref)) - ccall((:fftwf_flops,libfftwf), Void, - (PlanPtr,Ptr{Float64},Ptr{Float64},Ptr{Float64}), - plan, ptr, ptr + 8, ptr + 16) - (round(Int64, ref[][1]), round(Int64, ref[][2]), round(Int64, ref[][3])) -end -flops(plan::FFTWPlan) = let ops = arithmetic_ops(plan) - ops[1] + ops[2] + 2 * ops[3] # add + mul + 2*fma -end - -# Pretty-printing plans - -function showfftdims(io, sz::Dims, istride::Dims, T) - if isempty(sz) - print(io, "0-dimensional") - elseif length(sz) == 1 - print(io, sz[1], "-element") - else - print(io, join(sz, "×")) - end - if istride == colmajorstrides(sz) - print(io, " array of ", T) - else - print(io, " $istride-strided array of ", T) - end -end - -# The sprint_plan function was released in FFTW 3.3.4 -sprint_plan_(plan::FFTWPlan{<:fftwDouble}) = - ccall((:fftw_sprint_plan,libfftw), Ptr{UInt8}, (PlanPtr,), plan) -sprint_plan_(plan::FFTWPlan{<:fftwSingle}) = - ccall((:fftwf_sprint_plan,libfftwf), Ptr{UInt8}, (PlanPtr,), plan) -function sprint_plan(plan::FFTWPlan) - p = sprint_plan_(plan) - str = unsafe_string(p) - Libc.free(p) - return str -end - -function show(io::IO, p::cFFTWPlan{T,K,inplace}) where {T,K,inplace} - print(io, inplace ? "FFTW in-place " : "FFTW ", - K < 0 ? "forward" : "backward", " plan for ") - showfftdims(io, p.sz, p.istride, T) - version >= v"3.3.4" && print(io, "\n", sprint_plan(p)) -end - -function show(io::IO, p::rFFTWPlan{T,K,inplace}) where {T,K,inplace} - print(io, inplace ? "FFTW in-place " : "FFTW ", - K < 0 ? "real-to-complex" : "complex-to-real", - " plan for ") - showfftdims(io, p.sz, p.istride, T) - version >= v"3.3.4" && print(io, "\n", sprint_plan(p)) -end - -function show(io::IO, p::r2rFFTWPlan{T,K,inplace}) where {T,K,inplace} - print(io, inplace ? "FFTW in-place r2r " : "FFTW r2r ") - if isempty(K) - print(io, "0-dimensional") - elseif K == ntuple(i -> K[1], length(K)) - print(io, kind2string(K[1])) - if length(K) > 1 - print(io, "^", length(K)) - end - else - print(io, join(map(kind2string, K), "×")) - end - print(io, " plan for ") - showfftdims(io, p.sz, p.istride, T) - version >= v"3.3.4" && print(io, "\n", sprint_plan(p)) -end - -# Check whether a FFTWPlan is applicable to a given input array, and -# throw an informative error if not: -function assert_applicable(p::FFTWPlan{T}, X::StridedArray{T}) where T - if size(X) != p.sz - throw(ArgumentError("FFTW plan applied to wrong-size array")) - elseif strides(X) != p.istride - throw(ArgumentError("FFTW plan applied to wrong-strides array")) - elseif alignment_of(X) != p.ialign || p.flags & UNALIGNED != 0 - throw(ArgumentError("FFTW plan applied to array with wrong memory alignment")) - end -end - -function assert_applicable(p::FFTWPlan{T,K,inplace}, X::StridedArray{T}, Y::StridedArray) where {T,K,inplace} - assert_applicable(p, X) - if size(Y) != p.osz - throw(ArgumentError("FFTW plan applied to wrong-size output")) - elseif strides(Y) != p.ostride - throw(ArgumentError("FFTW plan applied to wrong-strides output")) - elseif alignment_of(Y) != p.oalign || p.flags & UNALIGNED != 0 - throw(ArgumentError("FFTW plan applied to output with wrong memory alignment")) - elseif inplace != (pointer(X) == pointer(Y)) - throw(ArgumentError(string("FFTW ", - inplace ? "in-place" : "out-of-place", - " plan applied to ", - inplace ? "out-of-place" : "in-place", - " data"))) - end -end - -# strides for a column-major (Julia-style) array of size == sz -colmajorstrides(sz) = isempty(sz) ? () : (1,cumprod(Int[sz[1:end-1]...])...) - -# Execute - -unsafe_execute!(plan::FFTWPlan{<:fftwDouble}) = - ccall((:fftw_execute,libfftw), Void, (PlanPtr,), plan) - -unsafe_execute!(plan::FFTWPlan{<:fftwSingle}) = - ccall((:fftwf_execute,libfftwf), Void, (PlanPtr,), plan) - -unsafe_execute!(plan::cFFTWPlan{T}, - X::StridedArray{T}, Y::StridedArray{T}) where {T<:fftwDouble} = - ccall((:fftw_execute_dft,libfftw), Void, - (PlanPtr,Ptr{T},Ptr{T}), plan, X, Y) - -unsafe_execute!(plan::cFFTWPlan{T}, - X::StridedArray{T}, Y::StridedArray{T}) where {T<:fftwSingle} = - ccall((:fftwf_execute_dft,libfftwf), Void, - (PlanPtr,Ptr{T},Ptr{T}), plan, X, Y) - -unsafe_execute!(plan::rFFTWPlan{Float64,FORWARD}, - X::StridedArray{Float64}, Y::StridedArray{Complex128}) = - ccall((:fftw_execute_dft_r2c,libfftw), Void, - (PlanPtr,Ptr{Float64},Ptr{Complex128}), plan, X, Y) - -unsafe_execute!(plan::rFFTWPlan{Float32,FORWARD}, - X::StridedArray{Float32}, Y::StridedArray{Complex64}) = - ccall((:fftwf_execute_dft_r2c,libfftwf), Void, - (PlanPtr,Ptr{Float32},Ptr{Complex64}), plan, X, Y) - -unsafe_execute!(plan::rFFTWPlan{Complex128,BACKWARD}, - X::StridedArray{Complex128}, Y::StridedArray{Float64}) = - ccall((:fftw_execute_dft_c2r,libfftw), Void, - (PlanPtr,Ptr{Complex128},Ptr{Float64}), plan, X, Y) - -unsafe_execute!(plan::rFFTWPlan{Complex64,BACKWARD}, - X::StridedArray{Complex64}, Y::StridedArray{Float32}) = - ccall((:fftwf_execute_dft_c2r,libfftwf), Void, - (PlanPtr,Ptr{Complex64},Ptr{Float32}), plan, X, Y) - -unsafe_execute!(plan::r2rFFTWPlan{T}, - X::StridedArray{T}, Y::StridedArray{T}) where {T<:fftwDouble} = - ccall((:fftw_execute_r2r,libfftw), Void, - (PlanPtr,Ptr{T},Ptr{T}), plan, X, Y) - -unsafe_execute!(plan::r2rFFTWPlan{T}, - X::StridedArray{T}, Y::StridedArray{T}) where {T<:fftwSingle} = - ccall((:fftwf_execute_r2r,libfftwf), Void, - (PlanPtr,Ptr{T},Ptr{T}), plan, X, Y) - -# NOTE ON GC (garbage collection): -# The FFTWPlan has a finalizer so that gc will destroy the plan, -# which is necessary for gc to work with plan_fft. However, -# even when we are creating a single-use FFTWPlan [e.g. for fftn(x)], -# we intentionally do NOT call destroy_plan explicitly, and instead -# wait for garbage collection. The reason is that, in the common -# case where the user calls fft(x) a second time soon afterwards, -# if destroy_plan has not yet been called then FFTW will internally -# re-use the table of trigonometric constants from the first plan. - -# Compute dims and howmany for FFTW guru planner -function dims_howmany(X::StridedArray, Y::StridedArray, - sz::Array{Int,1}, region) - reg = [region...] - if length(unique(reg)) < length(reg) - throw(ArgumentError("each dimension can be transformed at most once")) - end - ist = [strides(X)...] - ost = [strides(Y)...] - dims = [sz[reg] ist[reg] ost[reg]]' - oreg = [1:ndims(X);] - oreg[reg] = 0 - oreg = filter(d -> d > 0, oreg) - howmany = [sz[oreg] ist[oreg] ost[oreg]]' - return (dims, howmany) -end - -# check & convert kinds into int32 array with same length as region -function fix_kinds(region, kinds) - if length(kinds) != length(region) - if length(kinds) > length(region) - throw(ArgumentError("too many transform kinds")) - else - if isempty(kinds) - throw(ArgumentError("must supply a transform kind")) - end - k = Vector{Int32}(length(region)) - k[1:length(kinds)] = [kinds...] - k[length(kinds)+1:end] = kinds[end] - kinds = k - end - else - kinds = Int32[kinds...] - end - for i = 1:length(kinds) - if kinds[i] < 0 || kinds[i] > 10 - throw(ArgumentError("invalid transform kind")) - end - end - return kinds -end - -# low-level FFTWPlan creation (for internal use in FFTW module) - -for (Tr,Tc,fftw,lib) in ((:Float64,:Complex128,"fftw",libfftw), - (:Float32,:Complex64,"fftwf",libfftwf)) - @eval function (::Type{cFFTWPlan{$Tc,K,inplace,N}})(X::StridedArray{$Tc,N}, - Y::StridedArray{$Tc,N}, - region, flags::Integer, timelimit::Real) where {K,inplace,N} - direction = K - set_timelimit($Tr, timelimit) - R = isa(region, Tuple) ? region : copy(region) - dims, howmany = dims_howmany(X, Y, [size(X)...], R) - plan = ccall(($(string(fftw,"_plan_guru64_dft")),$lib), - PlanPtr, - (Int32, Ptr{Int}, Int32, Ptr{Int}, - Ptr{$Tc}, Ptr{$Tc}, Int32, UInt32), - size(dims,2), dims, size(howmany,2), howmany, - X, Y, direction, flags) - set_timelimit($Tr, NO_TIMELIMIT) - if plan == C_NULL - error("FFTW could not create plan") # shouldn't normally happen - end - return cFFTWPlan{$Tc,K,inplace,N}(plan, flags, R, X, Y) - end - - @eval function (::Type{rFFTWPlan{$Tr,$FORWARD,inplace,N}})(X::StridedArray{$Tr,N}, - Y::StridedArray{$Tc,N}, - region, flags::Integer, timelimit::Real) where {inplace,N} - R = isa(region, Tuple) ? region : copy(region) - region = circshift([region...],-1) # FFTW halves last dim - set_timelimit($Tr, timelimit) - dims, howmany = dims_howmany(X, Y, [size(X)...], region) - plan = ccall(($(string(fftw,"_plan_guru64_dft_r2c")),$lib), - PlanPtr, - (Int32, Ptr{Int}, Int32, Ptr{Int}, - Ptr{$Tr}, Ptr{$Tc}, UInt32), - size(dims,2), dims, size(howmany,2), howmany, - X, Y, flags) - set_timelimit($Tr, NO_TIMELIMIT) - if plan == C_NULL - error("FFTW could not create plan") # shouldn't normally happen - end - return rFFTWPlan{$Tr,$FORWARD,inplace,N}(plan, flags, R, X, Y) - end - - @eval function (::Type{rFFTWPlan{$Tc,$BACKWARD,inplace,N}})(X::StridedArray{$Tc,N}, - Y::StridedArray{$Tr,N}, - region, flags::Integer, timelimit::Real) where {inplace,N} - R = isa(region, Tuple) ? region : copy(region) - region = circshift([region...],-1) # FFTW halves last dim - set_timelimit($Tr, timelimit) - dims, howmany = dims_howmany(X, Y, [size(Y)...], region) - plan = ccall(($(string(fftw,"_plan_guru64_dft_c2r")),$lib), - PlanPtr, - (Int32, Ptr{Int}, Int32, Ptr{Int}, - Ptr{$Tc}, Ptr{$Tr}, UInt32), - size(dims,2), dims, size(howmany,2), howmany, - X, Y, flags) - set_timelimit($Tr, NO_TIMELIMIT) - if plan == C_NULL - error("FFTW could not create plan") # shouldn't normally happen - end - return rFFTWPlan{$Tc,$BACKWARD,inplace,N}(plan, flags, R, X, Y) - end - - @eval function (::Type{r2rFFTWPlan{$Tr,ANY,inplace,N}})(X::StridedArray{$Tr,N}, - Y::StridedArray{$Tr,N}, - region, kinds, flags::Integer, - timelimit::Real) where {inplace,N} - R = isa(region, Tuple) ? region : copy(region) - knd = fix_kinds(region, kinds) - set_timelimit($Tr, timelimit) - dims, howmany = dims_howmany(X, Y, [size(X)...], region) - plan = ccall(($(string(fftw,"_plan_guru64_r2r")),$lib), - PlanPtr, - (Int32, Ptr{Int}, Int32, Ptr{Int}, - Ptr{$Tr}, Ptr{$Tr}, Ptr{Int32}, UInt32), - size(dims,2), dims, size(howmany,2), howmany, - X, Y, knd, flags) - set_timelimit($Tr, NO_TIMELIMIT) - if plan == C_NULL - error("FFTW could not create plan") # shouldn't normally happen - end - r2rFFTWPlan{$Tr,(map(Int,knd)...),inplace,N}(plan, flags, R, X, Y) - end - - # support r2r transforms of complex = transforms of real & imag parts - @eval function (::Type{r2rFFTWPlan{$Tc,ANY,inplace,N}})(X::StridedArray{$Tc,N}, - Y::StridedArray{$Tc,N}, - region, kinds, flags::Integer, - timelimit::Real) where {inplace,N} - R = isa(region, Tuple) ? region : copy(region) - knd = fix_kinds(region, kinds) - set_timelimit($Tr, timelimit) - dims, howmany = dims_howmany(X, Y, [size(X)...], region) - dims[2:3, 1:size(dims,2)] *= 2 - howmany[2:3, 1:size(howmany,2)] *= 2 - howmany = [howmany [2,1,1]] # append loop over real/imag parts - plan = ccall(($(string(fftw,"_plan_guru64_r2r")),$lib), - PlanPtr, - (Int32, Ptr{Int}, Int32, Ptr{Int}, - Ptr{$Tc}, Ptr{$Tc}, Ptr{Int32}, UInt32), - size(dims,2), dims, size(howmany,2), howmany, - X, Y, knd, flags) - set_timelimit($Tr, NO_TIMELIMIT) - if plan == C_NULL - error("FFTW could not create plan") # shouldn't normally happen - end - r2rFFTWPlan{$Tc,(map(Int,knd)...),inplace,N}(plan, flags, R, X, Y) - end - -end - -# Convert arrays of numeric types to FFTW-supported packed complex-float types -# (FIXME: is there a way to use the Julia promotion rules more cleverly here?) -fftwcomplex(X::StridedArray{<:fftwComplex}) = X -fftwcomplex(X::AbstractArray{T}) where {T<:fftwReal} = - copy!(Array{typeof(complex(zero(T)))}(size(X)), X) -fftwcomplex(X::AbstractArray{<:Real}) = copy!(Array{Complex128}(size(X)),X) -fftwcomplex(X::AbstractArray{<:Complex}) = copy!(Array{Complex128}(size(X)), X) -fftwfloat(X::StridedArray{<:fftwReal}) = X -fftwfloat(X::AbstractArray{<:Real}) = copy!(Array{Float64}(size(X)), X) -fftwfloat(X::AbstractArray{<:Complex}) = fftwcomplex(X) - -for (f,direction) in ((:fft,FORWARD), (:bfft,BACKWARD)) - plan_f = Symbol("plan_",f) - plan_f! = Symbol("plan_",f,"!") - idirection = -direction - @eval begin - function $plan_f(X::StridedArray{T,N}, region; - flags::Integer=ESTIMATE, - timelimit::Real=NO_TIMELIMIT) where {T<:fftwComplex,N} - cFFTWPlan{T,$direction,false,N}(X, fakesimilar(flags, X, T), - region, flags, timelimit) - end - - function $plan_f!(X::StridedArray{T,N}, region; - flags::Integer=ESTIMATE, - timelimit::Real=NO_TIMELIMIT) where {T<:fftwComplex,N} - cFFTWPlan{T,$direction,true,N}(X, X, region, flags, timelimit) - end - $plan_f(X::StridedArray{<:fftwComplex}; kws...) = - $plan_f(X, 1:ndims(X); kws...) - $plan_f!(X::StridedArray{<:fftwComplex}; kws...) = - $plan_f!(X, 1:ndims(X); kws...) - - function plan_inv(p::cFFTWPlan{T,$direction,inplace,N}) where {T<:fftwComplex,N,inplace} - X = Array{T}(p.sz) - Y = inplace ? X : fakesimilar(p.flags, X, T) - ScaledPlan(cFFTWPlan{T,$idirection,inplace,N}(X, Y, p.region, - p.flags, NO_TIMELIMIT), - normalization(X, p.region)) - end - end -end - -function A_mul_B!(y::StridedArray{T}, p::cFFTWPlan{T}, x::StridedArray{T}) where T - assert_applicable(p, x, y) - unsafe_execute!(p, x, y) - return y -end - -function *(p::cFFTWPlan{T,K,false}, x::StridedArray{T,N}) where {T,K,N} - assert_applicable(p, x) - y = Array{T}(p.osz)::Array{T,N} - unsafe_execute!(p, x, y) - return y -end - -function *(p::cFFTWPlan{T,K,true}, x::StridedArray{T}) where {T,K} - assert_applicable(p, x) - unsafe_execute!(p, x, x) - return x -end - -# rfft/brfft and planned variants. No in-place version for now. - -for (Tr,Tc) in ((:Float32,:Complex64),(:Float64,:Complex128)) - # Note: use $FORWARD and $BACKWARD below because of issue #9775 - @eval begin - function plan_rfft(X::StridedArray{$Tr,N}, region; - flags::Integer=ESTIMATE, - timelimit::Real=NO_TIMELIMIT) where N - osize = rfft_output_size(X, region) - Y = flags&ESTIMATE != 0 ? FakeArray($Tc,osize...) : Array{$Tc}(osize...) - rFFTWPlan{$Tr,$FORWARD,false,N}(X, Y, region, flags, timelimit) - end - - function plan_brfft(X::StridedArray{$Tc,N}, d::Integer, region; - flags::Integer=ESTIMATE, - timelimit::Real=NO_TIMELIMIT) where N - osize = brfft_output_size(X, d, region) - Y = flags&ESTIMATE != 0 ? FakeArray($Tr,osize...) : Array{$Tr}(osize...) - - # FFTW currently doesn't support PRESERVE_INPUT for - # multidimensional out-of-place c2r transforms, so - # we have to handle 1d and >1d cases separately with a copy. Ugh. - if length(region) <= 1 - rFFTWPlan{$Tc,$BACKWARD,false,N}(X, Y, region, - flags | PRESERVE_INPUT, - timelimit) - else - rFFTWPlan{$Tc,$BACKWARD,false,N}(copy(X), Y, region, flags, - timelimit) - end - end - - plan_rfft(X::StridedArray{$Tr};kws...)=plan_rfft(X,1:ndims(X);kws...) - plan_brfft(X::StridedArray{$Tr};kws...)=plan_brfft(X,1:ndims(X);kws...) - - function plan_inv(p::rFFTWPlan{$Tr,$FORWARD,false,N}) where N - X = Array{$Tr}(p.sz) - Y = p.flags&ESTIMATE != 0 ? FakeArray($Tc,p.osz) : Array{$Tc}(p.osz) - ScaledPlan(rFFTWPlan{$Tc,$BACKWARD,false,N}(Y, X, p.region, - length(p.region) <= 1 ? - p.flags | PRESERVE_INPUT : - p.flags, NO_TIMELIMIT), - normalization(X, p.region)) - end - - function plan_inv(p::rFFTWPlan{$Tc,$BACKWARD,false,N}) where N - X = Arra{$Tc}(p.sz) - Y = p.flags&ESTIMATE != 0 ? FakeArray($Tr,p.osz) : Array{$Tr}(p.osz) - ScaledPlan(rFFTWPlan{$Tr,$FORWARD,false,N}(Y, X, p.region, - p.flags, NO_TIMELIMIT), - normalization(Y, p.region)) - end - - function A_mul_B!(y::StridedArray{$Tc}, p::rFFTWPlan{$Tr,$FORWARD}, x::StridedArray{$Tr}) - assert_applicable(p, x, y) - unsafe_execute!(p, x, y) - return y - end - function A_mul_B!(y::StridedArray{$Tr}, p::rFFTWPlan{$Tc,$BACKWARD}, x::StridedArray{$Tc}) - assert_applicable(p, x, y) - unsafe_execute!(p, x, y) # note: may overwrite x as well as y! - return y - end - - function *(p::rFFTWPlan{$Tr,$FORWARD,false}, x::StridedArray{$Tr,N}) where N - assert_applicable(p, x) - y = Array{$Tc}(p.osz)::Array{$Tc,N} - unsafe_execute!(p, x, y) - return y - end - - function *(p::rFFTWPlan{$Tc,$BACKWARD,false}, x::StridedArray{$Tc,N}) where N - if p.flags & PRESERVE_INPUT != 0 - assert_applicable(p, x) - y = Array{$Tr}(p.osz)::Array{$Tr,N} - unsafe_execute!(p, x, y) - else # need to make a copy to avoid overwriting x - xc = copy(x) - assert_applicable(p, xc) - y = Array{$Tr}(p.osz)::Array{$Tr,N} - unsafe_execute!(p, xc, y) - end - return y - end - end -end - -# FFTW r2r transforms (low-level interface) - -for f in (:r2r, :r2r!) - pf = Symbol("plan_", f) - @eval begin - $f(x::AbstractArray{<:fftwNumber}, kinds) = $pf(x, kinds) * x - $f(x::AbstractArray{<:fftwNumber}, kinds, region) = $pf(x, kinds, region) * x - $pf(x::AbstractArray, kinds; kws...) = $pf(x, kinds, 1:ndims(x); kws...) - $f(x::AbstractArray{<:Real}, kinds, region=1:ndims(x)) = $f(fftwfloat(x), kinds, region) - $pf(x::AbstractArray{<:Real}, kinds, region; kws...) = $pf(fftwfloat(x), kinds, region; kws...) - $f(x::AbstractArray{<:Complex}, kinds, region=1:ndims(x)) = $f(fftwcomplex(x), kinds, region) - $pf(x::AbstractArray{<:Complex}, kinds, region; kws...) = $pf(fftwcomplex(x), kinds, region; kws...) - end -end - -function plan_r2r(X::StridedArray{T,N}, kinds, region; - flags::Integer=ESTIMATE, - timelimit::Real=NO_TIMELIMIT) where {T<:fftwNumber,N} - r2rFFTWPlan{T,ANY,false,N}(X, fakesimilar(flags, X, T), region, kinds, - flags, timelimit) -end - -function plan_r2r!(X::StridedArray{T,N}, kinds, region; - flags::Integer=ESTIMATE, - timelimit::Real=NO_TIMELIMIT) where {T<:fftwNumber,N} - r2rFFTWPlan{T,ANY,true,N}(X, X, region, kinds, flags, timelimit) -end - -# mapping from r2r kind to the corresponding inverse transform -const inv_kind = Dict{Int,Int}(R2HC => HC2R, HC2R => R2HC, DHT => DHT, - REDFT00 => REDFT00, - REDFT01 => REDFT10, REDFT10 => REDFT01, - REDFT11 => REDFT11, - RODFT00 => RODFT00, - RODFT01 => RODFT10, RODFT10 => RODFT01, - RODFT11 => RODFT11) - -# r2r inverses are normalized to 1/N, where N is a "logical" size -# the transform with length n and kind k: -function logical_size(n::Integer, k::Integer) - k <= DHT && return n - k == REDFT00 && return 2(n-1) - k == RODFT00 && return 2(n+1) - return 2n -end - -function plan_inv(p::r2rFFTWPlan{T,K,inplace,N}) where {T<:fftwNumber,K,inplace,N} - X = Array{T}(p.sz) - iK = fix_kinds(p.region, [inv_kind[k] for k in K]) - Y = inplace ? X : fakesimilar(p.flags, X, T) - ScaledPlan(r2rFFTWPlan{T,ANY,inplace,N}(X, Y, p.region, iK, - p.flags, NO_TIMELIMIT), - normalization(real(T), - map(logical_size, [p.sz...][[p.region...]], iK), - 1:length(iK))) -end - -function A_mul_B!(y::StridedArray{T}, p::r2rFFTWPlan{T}, x::StridedArray{T}) where T - assert_applicable(p, x, y) - unsafe_execute!(p, x, y) - return y -end - -function *(p::r2rFFTWPlan{T,K,false}, x::StridedArray{T,N}) where {T,K,N} - assert_applicable(p, x) - y = Array{T}(p.osz)::Array{T,N} - unsafe_execute!(p, x, y) - return y -end - -function *(p::r2rFFTWPlan{T,K,true}, x::StridedArray{T}) where {T,K} - assert_applicable(p, x) - unsafe_execute!(p, x, x) - return x -end - -include("dct.jl") diff --git a/base/fft/dct.jl b/base/fft/dct.jl deleted file mode 100644 index a3410961a7f89..0000000000000 --- a/base/fft/dct.jl +++ /dev/null @@ -1,103 +0,0 @@ -# This file is a part of Julia. License is MIT: https://julialang.org/license - -# (This is part of the FFTW module.) - -export dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct! - -# Discrete cosine transforms (type II/III) via FFTW's r2r transforms; -# we follow the Matlab convention and adopt a unitary normalization here. -# Unlike Matlab we compute the multidimensional transform by default, -# similar to the Julia fft functions. - -mutable struct DCTPlan{T<:fftwNumber,K,inplace} <: Plan{T} - plan::r2rFFTWPlan{T} - r::Array{UnitRange{Int}} # array of indices for rescaling - nrm::Float64 # normalization factor - region::Dims # dimensions being transformed - pinv::DCTPlan{T} - DCTPlan{T,K,inplace}(plan,r,nrm,region) where {T<:fftwNumber,K,inplace} = new(plan,r,nrm,region) -end - -size(p::DCTPlan) = size(p.plan) - -function show(io::IO, p::DCTPlan{T,K,inplace}) where {T,K,inplace} - print(io, inplace ? "FFTW in-place " : "FFTW ", - K == REDFT10 ? "DCT (DCT-II)" : "IDCT (DCT-III)", " plan for ") - showfftdims(io, p.plan.sz, p.plan.istride, eltype(p)) -end - -for (pf, pfr, K, inplace) in ((:plan_dct, :plan_r2r, REDFT10, false), - (:plan_dct!, :plan_r2r!, REDFT10, true), - (:plan_idct, :plan_r2r, REDFT01, false), - (:plan_idct!, :plan_r2r!, REDFT01, true)) - @eval function $pf(X::StridedArray{T}, region; kws...) where T<:fftwNumber - r = [1:n for n in size(X)] - nrm = sqrt(0.5^length(region) * normalization(X,region)) - DCTPlan{T,$K,$inplace}($pfr(X, $K, region; kws...), r, nrm, - ntuple(i -> Int(region[i]), length(region))) - end -end - -function plan_inv(p::DCTPlan{T,K,inplace}) where {T,K,inplace} - X = Array{T}(p.plan.sz) - iK = inv_kind[K] - DCTPlan{T,iK,inplace}(inplace ? - plan_r2r!(X, iK, p.region, flags=p.plan.flags) : - plan_r2r(X, iK, p.region, flags=p.plan.flags), - p.r, p.nrm, p.region) -end - -for f in (:dct, :dct!, :idct, :idct!) - pf = Symbol("plan_", f) - @eval begin - $f(x::AbstractArray{<:fftwNumber}) = $pf(x) * x - $f(x::AbstractArray{<:fftwNumber}, region) = $pf(x, region) * x - $pf(x::AbstractArray; kws...) = $pf(x, 1:ndims(x); kws...) - $f(x::AbstractArray{<:Real}, region=1:ndims(x)) = $f(fftwfloat(x), region) - $pf(x::AbstractArray{<:Real}, region; kws...) = $pf(fftwfloat(x), region; kws...) - $pf(x::AbstractArray{<:Complex}, region; kws...) = $pf(fftwcomplex(x), region; kws...) - end -end - -const sqrthalf = sqrt(0.5) -const sqrt2 = sqrt(2.0) -const onerange = 1:1 - -function A_mul_B!(y::StridedArray{T}, p::DCTPlan{T,REDFT10}, - x::StridedArray{T}) where T - assert_applicable(p.plan, x, y) - unsafe_execute!(p.plan, x, y) - scale!(y, p.nrm) - r = p.r - for d in p.region - oldr = r[d] - r[d] = onerange - y[r...] *= sqrthalf - r[d] = oldr - end - return y -end - -# note: idct changes input data -function A_mul_B!(y::StridedArray{T}, p::DCTPlan{T,REDFT01}, - x::StridedArray{T}) where T - assert_applicable(p.plan, x, y) - scale!(x, p.nrm) - r = p.r - for d in p.region - oldr = r[d] - r[d] = onerange - x[r...] *= sqrt2 - r[d] = oldr - end - unsafe_execute!(p.plan, x, y) - return y -end - -*(p::DCTPlan{T,REDFT10,false}, x::StridedArray{T}) where {T} = - A_mul_B!(Array{T}(p.plan.osz), p, x) - -*(p::DCTPlan{T,REDFT01,false}, x::StridedArray{T}) where {T} = - A_mul_B!(Array{T}(p.plan.osz), p, copy(x)) # need copy to preserve input - -*(p::DCTPlan{T,K,true}, x::StridedArray{T}) where {T,K} = A_mul_B!(x, p, x) diff --git a/base/sysimg.jl b/base/sysimg.jl index acfb4c7b68ee4..a56d90a45d80e 100644 --- a/base/sysimg.jl +++ b/base/sysimg.jl @@ -337,12 +337,6 @@ include("statistics.jl") # irrational mathematical constants include("irrationals.jl") -# signal processing -include("dft.jl") -importall .DFT -include("dsp.jl") -importall .DSP - # Fast math include("fastmath.jl") importall .FastMath diff --git a/base/util.jl b/base/util.jl index b2167ca36b168..6b03ca0d8d48c 100644 --- a/base/util.jl +++ b/base/util.jl @@ -382,14 +382,6 @@ macro timed(ex) end end -function fftw_vendor() - if Base.libfftw_name in ("libmkl_rt", "mkl_rt") - return :mkl - else - return :fftw - end -end - ## printing with color ## diff --git a/contrib/julia.lang b/contrib/julia.lang index 80dff29ce6057..616e635de90af 100644 --- a/contrib/julia.lang +++ b/contrib/julia.lang @@ -219,7 +219,6 @@ Core Main - FFTW Collections Test Pkg diff --git a/contrib/julia.xml b/contrib/julia.xml index 978dc6e0b726e..577fc37de5d62 100644 --- a/contrib/julia.xml +++ b/contrib/julia.xml @@ -263,7 +263,6 @@ TypeError Collections - FFTW Order Sys Test diff --git a/contrib/mac/macports.make b/contrib/mac/macports.make index f99927ea90cfb..c4b32b2157d1b 100644 --- a/contrib/mac/macports.make +++ b/contrib/mac/macports.make @@ -22,7 +22,6 @@ JMAKEFLAGS = USE_SYSTEM_LLVM=1 LLVM_CONFIG=${LLVM_CONFIG} \ USE_SYSTEM_OPENLIBM=0 \ USE_SYSTEM_BLAS=1 USE_BLAS64=0\ USE_SYSTEM_LAPACK=1 \ - USE_SYSTEM_FFTW=1 \ USE_SYSTEM_GMP=1 \ USE_SYSTEM_MPFR=1 \ USE_SYSTEM_ARPACK=1 \ diff --git a/contrib/vagrant/Vagrantfile b/contrib/vagrant/Vagrantfile index 5ce593eb28655..d59ce3ee39a07 100644 --- a/contrib/vagrant/Vagrantfile +++ b/contrib/vagrant/Vagrantfile @@ -5,7 +5,7 @@ $script = <