Merge pull request #11 from stevengj/inplace

in-place natural-order transforms

README.md

@@ -51,7 +51,8 @@ and sums them to recover the signal, with no `n` factor.
   (Hadamard) ordering with `fwht_natural` and `ifwht_natural`, or in the
   dyadic (Paley) ordering with `fwht_dyadic` and `ifwht_dyadic`.  These
   functions take the same arguments as `fwht` and `ifwht` and have the
-  same normalizations, respectively.
+  same normalizations, respectively.    The natural-order transforms also
+  have in-place variants `fwht_natural!` and `ifwht_natural!`.
 
 ## Hadamard matrices
 

src/Hadamard.jl

@@ -15,7 +15,7 @@ There is also a function `hadamard(n)` that returns Hadamard matrices for known
 (including non powers of two).
 """
 module Hadamard
-export fwht, ifwht, fwht_natural, ifwht_natural, fwht_dyadic, ifwht_dyadic, hadamard
+export fwht, ifwht, fwht_natural, ifwht_natural, fwht_natural!, ifwht_natural!, fwht_dyadic, ifwht_dyadic, hadamard
 
 using FFTW
 import FFTW: set_timelimit, dims_howmany, unsafe_execute!, cFFTWPlan, r2rFFTWPlan, PlanPtr, fftwNumber, ESTIMATE, NO_TIMELIMIT, R2HC
@@ -205,6 +205,29 @@ function ifwht(X::StridedArray{<:Number}, region)
     return ifwht(float(X), region)
 end
 
+############################################################################
+# in-place transforms, currently for natural ordering only.
+
+"""
+    ifwht_natural!(X, dims=1:ndims(X))
+
+Similar to `ifwht_natural`, but works in-place on the input array `X`.
+"""
+function ifwht_natural!(X::StridedArray{<:fftwNumber}, region=1:ndims(X))
+    p = Plan_Hadamard(X, X, region, ESTIMATE, NO_TIMELIMIT, false)
+    unsafe_execute!(p)
+    return X
+end
+
+
+"""
+    fwht_natural!(X, dims=1:ndims(X))
+
+Similar to `fwht_natural`, but works in-place on the input array `X`.
+"""
+fwht_natural!(X::StridedArray{<:fftwNumber}, region=1:ndims(X)) =
+    Compat.rmul!(ifwht_natural!(X, region), normalization(X,region))
+
 ############################################################################
 # Forward transforms (normalized by 1/N as in Matlab) and transforms
 # without the region argument:

test/runtests.jl

@@ -109,6 +109,11 @@ end
 @test ifwht_natural(I32, 1) == hadamard(32)
 @test ifwht_natural(Matrix{Float64}(I,2,2), 1) == hadamard(2)
 
+let X = copy(I32)
+    @test ifwht_natural!(X, 1) == X == hadamard(32)
+    @test fwht_natural!(X, 1) == X == I32
+end
+
 print("Checking unitarity of hadamard(n): ")
 sizes = Int[]
 for i = 4:4:1200