From f286acf2ef8e76cd6f45fc88c51efcd57c51d3e4 Mon Sep 17 00:00:00 2001 From: Alex Arslan Date: Wed, 7 Jun 2017 16:23:58 -0700 Subject: [PATCH] Remove DSP functions (to be rehomed) (#11) * Remove DSP functions (to be rehomed) * Remove DSP doc file from the doc build * Copy necessary DCT tests to runtests.jl * Ensure DCT functions are called from the package, not Base --- docs/make.jl | 3 +- docs/src/dsp.md | 12 --- src/FFTW.jl | 4 - src/dsp.jl | 201 ----------------------------------------------- test/dsp.jl | 139 -------------------------------- test/runtests.jl | 95 ++++++++++++++++++++-- 6 files changed, 91 insertions(+), 363 deletions(-) delete mode 100644 docs/src/dsp.md delete mode 100644 src/dsp.jl delete mode 100644 test/dsp.jl diff --git a/docs/make.jl b/docs/make.jl index 31d9ce6..01ed471 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,8 +7,7 @@ makedocs( sitename = "FFTW.jl", pages = Any[ "Home" => "index.md", - "Fourier Transforms" => "fft.md", - "Signal Processing" => "dsp.md", + "API" => "fft.md", ], ) diff --git a/docs/src/dsp.md b/docs/src/dsp.md deleted file mode 100644 index 26e1bba..0000000 --- a/docs/src/dsp.md +++ /dev/null @@ -1,12 +0,0 @@ -# Signal Processing - -It is expected that these functions will be moved to a different package. - -```@docs -FFTW.filt -FFTW.filt! -FFTW.deconv -FFTW.conv -FFTW.conv2 -FFTW.xcorr -``` diff --git a/src/FFTW.jl b/src/FFTW.jl index 3ba8531..ec82703 100644 --- a/src/FFTW.jl +++ b/src/FFTW.jl @@ -16,9 +16,6 @@ import AbstractFFTs: Plan, ScaledPlan, if !isdefined(Base, :FFTW) export dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct! end -if !isdefined(Base, :DSP) - export filt, filt!, deconv, conv, conv2, xcorr -end const depsfile = joinpath(dirname(@__DIR__), "deps", "deps.jl") if isfile(depsfile) @@ -58,6 +55,5 @@ end include("fft.jl") include("dct.jl") -include("dsp.jl") # TODO: Move these functions to DSP.jl end # module diff --git a/src/dsp.jl b/src/dsp.jl deleted file mode 100644 index 34349d0..0000000 --- a/src/dsp.jl +++ /dev/null @@ -1,201 +0,0 @@ -# This file was formerly a part of Julia. License is MIT: https://julialang.org/license - -import Base.trailingsize - -_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 diff --git a/test/dsp.jl b/test/dsp.jl deleted file mode 100644 index b1af062..0000000 --- a/test/dsp.jl +++ /dev/null @@ -1,139 +0,0 @@ -# This file was formerly a part of Julia. License is MIT: https://julialang.org/license - -# Filter -b = [1., 2., 3., 4.] -x = [1., 1., 0., 1., 1., 0., 0., 0.] -@test filt(b, 1., x) == [1., 3., 5., 8., 7., 5., 7., 4.] -@test filt(b, [1., -0.5], x) == [1., 3.5, 6.75, 11.375, 12.6875, 11.34375, 12.671875, 10.3359375] -# With ranges -@test filt(b, 1., 1.0:10.0) == [1., 4., 10., 20., 30., 40., 50., 60., 70., 80.] -@test filt(1.:4., 1., 1.0:10.0) == [1., 4., 10., 20., 30., 40., 50., 60., 70., 80.] -# Across an array is the same as channel-by-channel -@test filt(b, 1., [x 1.0:8.0]) == [filt(b, 1., x) filt(b, 1., 1.0:8.0)] -@test filt(b, [1., -0.5], [x 1.0:8.0]) == [filt(b, [1., -0.5], x) filt(b, [1., -0.5], 1.0:8.0)] -si = zeros(3) -@test filt(b, 1., [x 1.0:8.0], si) == [filt(b, 1., x, si) filt(b, 1., 1.0:8.0, si)] -@test si == zeros(3) # Will likely fail if/when arrayviews are implemented -si = [zeros(3) ones(3)] -@test filt(b, 1., [x 1.0:8.0], si) == [filt(b, 1., x, zeros(3)) filt(b, 1., 1.0:8.0, ones(3))] -# With initial conditions: a lowpass 5-pole butterworth filter with W_n = 0.25, -# and a stable initial filter condition matched to the initial value -b = [0.003279216306360201,0.016396081531801006,0.03279216306360201,0.03279216306360201,0.016396081531801006,0.003279216306360201] -a = [1.0,-2.4744161749781606,2.8110063119115782,-1.703772240915465,0.5444326948885326,-0.07231566910295834] -si = [0.9967207836936347,-1.4940914728163142,1.2841226760316475,-0.4524417279474106,0.07559488540931815] -@test filt(b, a, ones(10), si) ≈ ones(10) # Shouldn't affect DC offset - -@test_throws ArgumentError filt!([1, 2], [1], [1], [1]) -@test xcorr([1, 2], [3, 4]) == [4, 11, 6] - -# Shift-Functions -@test fftshift([1 2 3]) == [3 1 2] -@test fftshift([1, 2, 3]) == [3, 1, 2] -@test fftshift([1 2 3; 4 5 6]) == [6 4 5; 3 1 2] - -@test fftshift([1 2 3; 4 5 6], 1) == [4 5 6; 1 2 3] -@test fftshift([1 2 3; 4 5 6], ()) == [1 2 3; 4 5 6] -@test fftshift([1 2 3; 4 5 6], (1,2)) == [6 4 5; 3 1 2] -@test fftshift([1 2 3; 4 5 6], 1:2) == [6 4 5; 3 1 2] - -@test ifftshift([1 2 3]) == [2 3 1] -@test ifftshift([1, 2, 3]) == [2, 3, 1] -@test ifftshift([1 2 3; 4 5 6]) == [5 6 4; 2 3 1] - -@test ifftshift([1 2 3; 4 5 6], 1) == [4 5 6; 1 2 3] -@test ifftshift([1 2 3; 4 5 6], ()) == [1 2 3; 4 5 6] -@test ifftshift([1 2 3; 4 5 6], (1,2)) == [5 6 4; 2 3 1] -@test ifftshift([1 2 3; 4 5 6], 1:2) == [5 6 4; 2 3 1] - -# Convolution -a = [1., 2., 1., 2.] -b = [1., 2., 3.] -@test conv(a, b) ≈ [1., 4., 8., 10., 7., 6.] -@test conv(complex.(a, ones(4)), complex(b)) ≈ complex.([1., 4., 8., 10., 7., 6.], [1., 3., 6., 6., 5., 3.]) - -# Discrete cosine transform (DCT) tests - -if FFTW.fftw_vendor() != :mkl - a = rand(8,11) + im*rand(8,11) - @test norm(idct(dct(a)) - a) < 1e-8 - - X = reshape([1,2,7,2,1,5,9,-1,3,4,6,9],3,4) - Y = rand(17,14) - Y[3:5,9:12] = X - sX = view(Y,3:5,9:12) - - true_Xdct = [ 13.856406460551018 -3.863239728836245 2.886751345948129 -0.274551994240164; -2.828427124746190 -2.184015211898548 -4.949747468305834 3.966116180118245; 4.898979485566356 -0.194137576915510 -2.857738033247041 2.731723009609389 ] - - true_Xdct_1 = [ 5.773502691896258 4.618802153517007 6.350852961085884 10.969655114602890; -4.242640687119286 -2.121320343559643 4.242640687119286 -3.535533905932738; 1.632993161855452 2.041241452319315 5.715476066494083 0.408248290463863 ] - - true_Xdct_2 = [ 8. -3.854030797826254 -3.0 3.761176226848022; - 4.0 -2.071929829606556 4.0 -2.388955165168770; 12. -0.765366864730179 4.0 -1.847759065022573 ] - - Xdct = dct(X) - Xdct! = float(X); dct!(Xdct!) - Xdct_1 = dct(X,1) - Xdct!_1 = float(X); dct!(Xdct!_1,1) - Xdct_2 = dct(X,2) - Xdct!_2 = float(X); dct!(Xdct!_2,2) - - Xidct = idct(true_Xdct) - Xidct! = copy(true_Xdct); idct!(Xidct!) - Xidct_1 = idct(true_Xdct_1,1) - Xidct!_1 = copy(true_Xdct_1); idct!(Xidct!_1,1) - Xidct_2 = idct(true_Xdct_2,2) - Xidct!_2 = copy(true_Xdct_2); idct!(Xidct!_2,2) - - pXdct = plan_dct(X)*(X) - pXdct! = float(X); plan_dct!(pXdct!)*(pXdct!) - pXdct_1 = plan_dct(X,1)*(X) - pXdct!_1 = float(X); plan_dct!(pXdct!_1,1)*(pXdct!_1) - pXdct_2 = plan_dct(X,2)*(X) - pXdct!_2 = float(X); plan_dct!(pXdct!_2,2)*(pXdct!_2) - - pXidct = plan_idct(true_Xdct)*(true_Xdct) - pXidct! = copy(true_Xdct); plan_idct!(pXidct!)*(pXidct!) - pXidct_1 = plan_idct(true_Xdct_1,1)*(true_Xdct_1) - pXidct!_1 = copy(true_Xdct_1); plan_idct!(pXidct!_1,1)*(pXidct!_1) - pXidct_2 = plan_idct(true_Xdct_2,2)*(true_Xdct_2) - pXidct!_2 = copy(true_Xdct_2); plan_idct!(pXidct!_2,2)*(pXidct!_2) - - sXdct = dct(sX) - psXdct = plan_dct(sX)*(sX) - sYdct! = copy(Y); sXdct! = view(sYdct!,3:5,9:12); dct!(sXdct!) - psYdct! = copy(Y); psXdct! = view(psYdct!,3:5,9:12); plan_dct!(psXdct!)*(psXdct!) - - for i = 1:length(X) - @test Xdct[i] ≈ true_Xdct[i] - @test Xdct![i] ≈ true_Xdct[i] - @test Xdct_1[i] ≈ true_Xdct_1[i] - @test Xdct!_1[i] ≈ true_Xdct_1[i] - @test Xdct_2[i] ≈ true_Xdct_2[i] - @test Xdct!_2[i] ≈ true_Xdct_2[i] - - @test pXdct[i] ≈ true_Xdct[i] - @test pXdct![i] ≈ true_Xdct[i] - @test pXdct_1[i] ≈ true_Xdct_1[i] - @test pXdct!_1[i] ≈ true_Xdct_1[i] - @test pXdct_2[i] ≈ true_Xdct_2[i] - @test pXdct!_2[i] ≈ true_Xdct_2[i] - - @test Xidct[i] ≈ X[i] - @test Xidct![i] ≈ X[i] - @test Xidct_1[i] ≈ X[i] - @test Xidct!_1[i] ≈ X[i] - @test Xidct_2[i] ≈ X[i] - @test Xidct!_2[i] ≈ X[i] - - @test pXidct[i] ≈ X[i] - @test pXidct![i] ≈ X[i] - @test pXidct_1[i] ≈ X[i] - @test pXidct!_1[i] ≈ X[i] - @test pXidct_2[i] ≈ X[i] - @test pXidct!_2[i] ≈ X[i] - - @test sXdct[i] ≈ true_Xdct[i] - @test psXdct[i] ≈ true_Xdct[i] - @test sXdct![i] ≈ true_Xdct[i] - @test psXdct![i] ≈ true_Xdct[i] - end -end # fftw_vendor() != :mkl diff --git a/test/runtests.jl b/test/runtests.jl index d24def4..dbd309f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,8 +8,9 @@ import AbstractFFTs: Plan, 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, fftshift, ifftshift, plan_inv +import FFTW: dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct! -# issue #19892 +# Base Julia issue #19892 # (test this first to make sure it happens before set_num_threads) let a = randn(10^5,1), p1 = plan_rfft(a, flags=FFTW.ESTIMATE) FFTW.set_num_threads(2) @@ -326,7 +327,7 @@ let FFTW.flops(plan64) end -# issue #9772 +# Base Julia issue #9772 for x in (randn(10),randn(10,12)) z = complex(x) y = rfft(x) @@ -352,7 +353,7 @@ for x in (randn(10),randn(10,12)) # algorithm steps are included in the CTPlan type end -# issue #17896 +# Base Julia issue #17896 a = rand(5) @test fft(a) == fft(view(a,:)) == fft(view(a, 1:5)) == fft(view(a, [1:5;])) @test rfft(a) == rfft(view(a,:)) == rfft(view(a, 1:5)) == rfft(view(a, [1:5;])) @@ -360,5 +361,89 @@ a16 = convert(Vector{Float16}, a) @test fft(a16) == fft(view(a16,:)) == fft(view(a16, 1:5)) == fft(view(a16, [1:5;])) @test rfft(a16) == rfft(view(a16,:)) == rfft(view(a16, 1:5)) == rfft(view(a16, [1:5;])) -# TODO: Move DSP functionality to DSP.jl -include("dsp.jl") +# Discrete cosine transform (DCT) tests + +if FFTW.fftw_vendor() != :mkl + a = rand(8,11) + im*rand(8,11) + @test norm(idct(dct(a)) - a) < 1e-8 + + X = reshape([1,2,7,2,1,5,9,-1,3,4,6,9],3,4) + Y = rand(17,14) + Y[3:5,9:12] = X + sX = view(Y,3:5,9:12) + + true_Xdct = [ 13.856406460551018 -3.863239728836245 2.886751345948129 -0.274551994240164; -2.828427124746190 -2.184015211898548 -4.949747468305834 3.966116180118245; 4.898979485566356 -0.194137576915510 -2.857738033247041 2.731723009609389 ] + + true_Xdct_1 = [ 5.773502691896258 4.618802153517007 6.350852961085884 10.969655114602890; -4.242640687119286 -2.121320343559643 4.242640687119286 -3.535533905932738; 1.632993161855452 2.041241452319315 5.715476066494083 0.408248290463863 ] + + true_Xdct_2 = [ 8. -3.854030797826254 -3.0 3.761176226848022; + 4.0 -2.071929829606556 4.0 -2.388955165168770; 12. -0.765366864730179 4.0 -1.847759065022573 ] + + Xdct = dct(X) + Xdct! = float(X); dct!(Xdct!) + Xdct_1 = dct(X,1) + Xdct!_1 = float(X); dct!(Xdct!_1,1) + Xdct_2 = dct(X,2) + Xdct!_2 = float(X); dct!(Xdct!_2,2) + + Xidct = idct(true_Xdct) + Xidct! = copy(true_Xdct); idct!(Xidct!) + Xidct_1 = idct(true_Xdct_1,1) + Xidct!_1 = copy(true_Xdct_1); idct!(Xidct!_1,1) + Xidct_2 = idct(true_Xdct_2,2) + Xidct!_2 = copy(true_Xdct_2); idct!(Xidct!_2,2) + + pXdct = plan_dct(X)*(X) + pXdct! = float(X); plan_dct!(pXdct!)*(pXdct!) + pXdct_1 = plan_dct(X,1)*(X) + pXdct!_1 = float(X); plan_dct!(pXdct!_1,1)*(pXdct!_1) + pXdct_2 = plan_dct(X,2)*(X) + pXdct!_2 = float(X); plan_dct!(pXdct!_2,2)*(pXdct!_2) + + pXidct = plan_idct(true_Xdct)*(true_Xdct) + pXidct! = copy(true_Xdct); plan_idct!(pXidct!)*(pXidct!) + pXidct_1 = plan_idct(true_Xdct_1,1)*(true_Xdct_1) + pXidct!_1 = copy(true_Xdct_1); plan_idct!(pXidct!_1,1)*(pXidct!_1) + pXidct_2 = plan_idct(true_Xdct_2,2)*(true_Xdct_2) + pXidct!_2 = copy(true_Xdct_2); plan_idct!(pXidct!_2,2)*(pXidct!_2) + + sXdct = dct(sX) + psXdct = plan_dct(sX)*(sX) + sYdct! = copy(Y); sXdct! = view(sYdct!,3:5,9:12); dct!(sXdct!) + psYdct! = copy(Y); psXdct! = view(psYdct!,3:5,9:12); plan_dct!(psXdct!)*(psXdct!) + + for i = 1:length(X) + @test Xdct[i] ≈ true_Xdct[i] + @test Xdct![i] ≈ true_Xdct[i] + @test Xdct_1[i] ≈ true_Xdct_1[i] + @test Xdct!_1[i] ≈ true_Xdct_1[i] + @test Xdct_2[i] ≈ true_Xdct_2[i] + @test Xdct!_2[i] ≈ true_Xdct_2[i] + + @test pXdct[i] ≈ true_Xdct[i] + @test pXdct![i] ≈ true_Xdct[i] + @test pXdct_1[i] ≈ true_Xdct_1[i] + @test pXdct!_1[i] ≈ true_Xdct_1[i] + @test pXdct_2[i] ≈ true_Xdct_2[i] + @test pXdct!_2[i] ≈ true_Xdct_2[i] + + @test Xidct[i] ≈ X[i] + @test Xidct![i] ≈ X[i] + @test Xidct_1[i] ≈ X[i] + @test Xidct!_1[i] ≈ X[i] + @test Xidct_2[i] ≈ X[i] + @test Xidct!_2[i] ≈ X[i] + + @test pXidct[i] ≈ X[i] + @test pXidct![i] ≈ X[i] + @test pXidct_1[i] ≈ X[i] + @test pXidct!_1[i] ≈ X[i] + @test pXidct_2[i] ≈ X[i] + @test pXidct!_2[i] ≈ X[i] + + @test sXdct[i] ≈ true_Xdct[i] + @test psXdct[i] ≈ true_Xdct[i] + @test sXdct![i] ≈ true_Xdct[i] + @test psXdct![i] ≈ true_Xdct[i] + end +end # fftw_vendor() != :mkl