Add walsh(n) function that returns a Walsh matrix (#28)

* Add walsh(n) function that returns a Walsh matrix by bit-reversal permutation of a Hadamard matrix.

* Apply suggestions from code review

Co-authored-by: Steven G. Johnson <[email protected]>

* Add tests for walsh(n)

* Update README.md for addition of walsth(n)

* Update README.md for addition of walsh(n)

---------

Co-authored-by: Steven G. Johnson <[email protected]>

Author: syoshida1983

README.md

@@ -68,6 +68,11 @@ choice is equivalent to `ifwht_natural(eye(n), 1)`.
 You can pretty-print a Hadamard matrix as a table of `+` and `-`
 (characters indicating the signs of the entries) via `Hadamard.printsigns`, e.g. `Hadamard.printsigns(hadamard(28))` for the 28×28 Hadamard matrix.
 
+You can also obtain a Walsh matrix (sequency-ordered Hadamard matrix) of
+order `n` by using the function `walsh(n)`; the order `n` *must* be powers
+of two. This function is related to the Hadamard matrix `hadamard(n)` by
+a bit-reversal permutation followed by a Gray-code permutation of the rows.
+
 ## Author
 
 This package was written by [Steven G. Johnson](http://math.mit.edu/~stevenj/).

src/Hadamard.jl

@@ -14,7 +14,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_natural!, ifwht_natural!, fwht_dyadic, ifwht_dyadic, hadamard
+export fwht, ifwht, fwht_natural, ifwht_natural, fwht_natural!, ifwht_natural!, fwht_dyadic, ifwht_dyadic, hadamard, walsh
 
 using FFTW, LinearAlgebra
 import FFTW: set_timelimit, dims_howmany, unsafe_execute!, cFFTWPlan, r2rFFTWPlan, PlanPtr, fftwNumber, ESTIMATE, NO_TIMELIMIT, R2HC
@@ -368,6 +368,33 @@ function hadamard(n::Integer)
     return H
 end
 
+"""
+    walsh(n)
+
+Return a Walsh matrix of order `n`, which must be a power of two, in sequency ordering. 
+This is related to the Hadamard matrix [`hadamard(n)`](@ref) by a bit-reversal permutation
+followed by a Gray-code permutation of the rows.
+
+Note that the linear operation `walsh(n) * x` can be computed much more efficiently, albeit
+in floating-point arithmetic, by `fwht(x) * n` (where the `* n` is due to the differing normalization,
+and can usually be eliminated by adjusting how you use the resulting vector) using the
+sequency-ordered fast Walsh–Hadamard transform function [`fwht`](@ref).
+"""
+walsh(n::Integer) = walsh(Int(n))
+
+function walsh(n::Int)
+    ispow2(n) || throw(ArgumentError("n=$n is not a power of two"))
+    m = trailing_zeros(n) + 1 # number of bits to represent the index
+    e = sizeof(Int)*8 - m       # number of extra bits
+    j = [ let b = bitreverse(i) >> e # bit-reverse the binary index (trailing m bits)
+              j = b ⊻ (b >> 1) # binary sequency index
+              j + 1 # 1-based index
+          end
+          for i in 0:n-1 ] 
+
+    return hadamard(n)[j, :]
+end
+
 ############################################################################
 
 end # module

test/runtests.jl

@@ -128,4 +128,13 @@ for i = 4:4:1200
     end
 end
 @test sizes == [4, 8, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48, 52, 56, 60, 64, 68, 72, 76, 80, 84, 88, 92, 96, 100, 104, 108, 112, 116, 120, 124, 128, 132, 136, 140, 144, 148, 152, 156, 160, 164, 168, 172, 176, 180, 184, 188, 192, 196, 200, 204, 208, 212, 216, 220, 224, 228, 232, 236, 240, 244, 248, 252, 256, 264, 272, 280, 288, 296, 304, 312, 320, 328, 336, 344, 352, 360, 368, 376, 384, 392, 400, 408, 416, 424, 428, 432, 440, 448, 456, 464, 472, 480, 488, 496, 504, 512, 528, 544, 560, 576, 592, 608, 624, 640, 656, 672, 688, 704, 720, 736, 752, 768, 784, 800, 816, 832, 848, 856, 864, 880, 896, 912, 928, 944, 960, 976, 992, 1008, 1024, 1040, 1056, 1088, 1104, 1120, 1152, 1184, 1200]
+
+@testset "walsh(n)" begin
+    @test walsh(8) == [1 1 1 1 1 1 1 1; 1 1 1 1 -1 -1 -1 -1; 1 1 -1 -1 -1 -1 1 1; 1 1 -1 -1 1 1 -1 -1; 1 -1 -1 1 1 -1 -1 1; 1 -1 -1 1 -1 1 1 -1; 1 -1 1 -1 -1 1 -1 1; 1 -1 1 -1 1 -1 1 -1]
+    @test_throws ArgumentError walsh(24)
+    for n in 2 .^ (0:10)
+        @test fwht(Matrix(I(n)), 1) * n ≈ walsh(n) rtol=1e-13
+    end
+end
+
 println(".\nSuccess!")