Skip to content

Commit cf97709

Browse files
stevengjstevengj
committed
initial working version
1 parent d5953f0 commit cf97709

File tree

4 files changed

+383
-0
lines changed

4 files changed

+383
-0
lines changed

‎LICENSE.md

+21
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
Copyright © 2013 by Steven G. Johnson
2+
3+
Permission is hereby granted, free of charge, to any person obtaining
4+
a copy of this software and associated documentation files (the
5+
"Software"), to deal in the Software without restriction, including
6+
without limitation the rights to use, copy, modify, merge, publish,
7+
distribute, sublicense, and/or sell copies of the Software, and to
8+
permit persons to whom the Software is furnished to do so, subject to
9+
the following conditions:
10+
11+
The above copyright notice and this permission notice shall be
12+
included in all copies or substantial portions of the Software.
13+
14+
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
15+
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
16+
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
17+
NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
18+
LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
19+
OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
20+
WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
21+

‎README.md

+58
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
# Fast multidimensional Walsh-Hadamard transforms
2+
3+
This package provides functions to compute fast Walsh-Hadamard transforms
4+
in Julia, for arbitrary dimensions and arbitrary power-of-two transform sizes,
5+
with the three standard orderings: natural (Hadamard), dyadic (Paley), and
6+
sequency (Walsh) ordering.
7+
8+
It works by calling Julia's interface to the [FFTW](http://www.fftw.org/)
9+
library, and can often be orders of magnitude faster than the corresponding
10+
`fwht` functions in the Matlab signal-processing toolbox.
11+
12+
## Installation
13+
14+
Within Julia, just use the package manager to run `Pkg.add("Hadamard")` to
15+
install the files.
16+
17+
## Usage
18+
19+
After installation, the `using Hadamard` statement will import the names
20+
in the Hadamard module so that you can call the function below.
21+
22+
* The function `fwht(X)` computes the fast Walsh-Hadamard transform
23+
(WHT) of the multidimensional array `X` (of real or complex numbers),
24+
returning its output in sequency order. The inverse transform is
25+
`ifwht(X)`.
26+
27+
By default, `fwht` and `ifwht` compute the *multidimensional* WHT, which
28+
consists of the ordinary (one-dimensional) WHT performed along each dimension
29+
of the input. To perform only the 1d WHT along dimension `d`, you can
30+
can instead use `fwht(X, d)` and `ifwht(X, d)` functions. More generally,
31+
`d` can be a tuple or array or dimensions to transform.
32+
33+
The sizes of the transformed dimensions *must* be powers of two, or an
34+
exception is thrown. The non-transformed dimensions are arbitrary. For
35+
example, given a 16x20 array `X`, `fwht(X,1)` is allowed but `fwht(X,2)` is
36+
not.
37+
38+
These functions compute the WHT normalized similarly to the `fwht` and
39+
`ifwht` functions in Matlab. Given the Walsh functions, which have values
40+
of +1 or -1, `fwht` multiplies its input by the Walsh functions and divides
41+
by `n` (the length of the input) to obtain the coefficients of each Walsh
42+
function in the input. `ifwht` multiplies its inputs by the Walsh functions
43+
and sums them to recover the signal, with no `n` factor.
44+
45+
* Instead of sequency order, one can also compute the WHT in the natural
46+
(Hadamard) ordering with `fwht_natural` and `ifwht_natural`, or in the
47+
dyadic (Paley) ordering with `fwht_dyadic` and `ifwht_dyadic`. These
48+
functions take the same arguments as `fwht` and `ifwht` and have the
49+
same normalizations, respectively.
50+
51+
* Also provided is a function `hadamard(n)` which returns the
52+
(natural-order) Hadamard matrix of order `n`, similar to the Matlab
53+
function of the same name. Currently, `n` must be a power of two,
54+
in which case this function is equivalent to `ifwht_natural(eye(n), 1)`.
55+
56+
## Author
57+
58+
This package was written by [Steven G. Johnson](http://math.mit.edu/~stevenj/).

‎src/Hadamard.jl

+219
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,219 @@
1+
# Compute fast Walsh-Hadamard transforms (in natural, dyadic, or sequency order)
2+
# using FFTW, by interpreting them as 2x2x2...x2x2 DFTs. We follow Matlab's
3+
# convention in which ifwht is the unnormalized transform and fwht has a 1/N
4+
# normalization (as opposed to using a unitary normalization).
5+
6+
module Hadamard
7+
export fwht, ifwht, fwht_natural, ifwht_natural, fwht_dyadic, ifwht_dyadic, hadamard
8+
9+
importall Base.FFTW
10+
import FFTW.set_timelimit
11+
import FFTW.dims_howmany
12+
import FFTW.Plan
13+
import FFTW.execute
14+
import FFTW.complexfloat
15+
import FFTW.normalization
16+
17+
power_of_two(n::Integer) = n > 0 && (n & (n - 1)) == 0
18+
19+
# A power-of-two dimension to be transformed is interpreted as a
20+
# 2x2x2x....x2x2 multidimensional DFT. This function transforms
21+
# each individual power-of-two dimension n into the corresponding log2(n)
22+
# DFT dimensions. If bitreverse is true, then the output is in
23+
# bit-reversed (transposed) order (which may not work in-place).
24+
function hadamardize(dims::Array{Int,2}, bitreverse::Bool)
25+
ntot = 0
26+
for i = 1:size(dims,2)
27+
n = dims[1,i]
28+
if !power_of_two(n)
29+
throw(ArgumentError("non power-of-two Hadamard-transform length"))
30+
end
31+
ntot += trailing_zeros(n)
32+
end
33+
hdims = Array(Int,3,ntot)
34+
j = 0
35+
for i = 1:size(dims,2)
36+
n = dims[1,i]
37+
is = dims[2,i]
38+
os = dims[3,i]
39+
log2n = trailing_zeros(n)
40+
krange = j+1:j+log2n
41+
for k in krange
42+
hdims[1,k] = 2
43+
hdims[2,k] = is
44+
hdims[3,k] = os
45+
is *= 2
46+
os *= 2
47+
end
48+
if bitreverse
49+
hdims[3,krange] = fliplr(hdims[3,krange])
50+
end
51+
j += log2n
52+
end
53+
return hdims
54+
end
55+
56+
for (Tr,Tc,fftw,lib) in ((:Float64,:Complex128,"fftw",FFTW.libfftw),
57+
(:Float32,:Complex64,"fftwf",FFTW.libfftwf))
58+
@eval function Plan_Hadamard(X::StridedArray{$Tc}, Y::StridedArray{$Tc},
59+
region, flags::Unsigned, timelimit::Real,
60+
bitreverse::Bool)
61+
set_timelimit($Tr, timelimit)
62+
dims, howmany = dims_howmany(X, Y, [size(X)...], region)
63+
dims = hadamardize(dims, bitreverse)
64+
plan = ccall(($(string(fftw,"_plan_guru64_dft")),$lib),
65+
Ptr{Void},
66+
(Int32, Ptr{Int}, Int32, Ptr{Int},
67+
Ptr{$Tc}, Ptr{$Tc}, Int32, Uint32),
68+
size(dims,2), dims, size(howmany,2), howmany,
69+
X, Y, FFTW.FORWARD, flags)
70+
set_timelimit($Tr, NO_TIMELIMIT)
71+
if plan == C_NULL
72+
error("FFTW could not create plan") # shouldn't normally happen
73+
end
74+
return Plan(plan, X)
75+
end
76+
77+
@eval function Plan_Hadamard(X::StridedArray{$Tr}, Y::StridedArray{$Tr},
78+
region, flags::Unsigned, timelimit::Real,
79+
bitreverse::Bool)
80+
set_timelimit($Tr, timelimit)
81+
dims, howmany = dims_howmany(X, Y, [size(X)...], region)
82+
dims = hadamardize(dims, bitreverse)
83+
kind = Array(Int32, size(dims,2))
84+
kind[:] = R2HC
85+
plan = ccall(($(string(fftw,"_plan_guru64_r2r")),$lib),
86+
Ptr{Void},
87+
(Int32, Ptr{Int}, Int32, Ptr{Int},
88+
Ptr{$Tr}, Ptr{$Tr}, Ptr{Int32}, Uint32),
89+
size(dims,2), dims, size(howmany,2), howmany,
90+
X, Y, kind, flags)
91+
set_timelimit($Tr, NO_TIMELIMIT)
92+
if plan == C_NULL
93+
error("FFTW could not create plan") # shouldn't normally happen
94+
end
95+
return Plan(plan, X)
96+
end
97+
end
98+
99+
############################################################################
100+
# Define ifwht (inverse/unnormalized) transforms for various orderings
101+
102+
# Natural (Hadamard) ordering:
103+
104+
function ifwht_natural{T<:fftwNumber}(X::StridedArray{T}, region)
105+
Y = similar(X)
106+
p = Plan_Hadamard(X, Y, region, ESTIMATE, NO_TIMELIMIT, false)
107+
execute(T, p.plan)
108+
return Y
109+
end
110+
111+
function ifwht_natural{T<:Number}(X::StridedArray{T}, region)
112+
Y = T<:Complex ? complexfloat(X) : float(X)
113+
p = Plan_Hadamard(Y, Y, region, ESTIMATE, NO_TIMELIMIT, false)
114+
execute(p.plan, Y, Y)
115+
return Y
116+
end
117+
118+
# Dyadic (Paley, bit-reversed) ordering:
119+
120+
function ifwht_dyadic{T<:fftwNumber}(X::StridedArray{T}, region)
121+
Y = similar(X)
122+
p = Plan_Hadamard(X, Y, region, ESTIMATE, NO_TIMELIMIT, true)
123+
execute(T, p.plan)
124+
return Y
125+
end
126+
127+
function ifwht_dyadic{T<:Number}(X::StridedArray{T}, region)
128+
X = T<:Complex ? complexfloat(X) : float(X)
129+
Y = similar(X) # FFTW cannot do bit-reversal in-place
130+
p = Plan_Hadamard(X, Y, region, ESTIMATE, NO_TIMELIMIT, true)
131+
execute(p.plan, X, Y)
132+
return Y
133+
end
134+
135+
############################################################################
136+
# Sequency (Walsh) ordering:
137+
138+
# ifwht along a single dimension d of X
139+
function ifwht{T<:fftwNumber}(X::Array{T}, region)
140+
Y = ifwht_dyadic(X, region)
141+
142+
# Perform Gray-code permutation of Y (TODO: in-place?)
143+
if isempty(region)
144+
return Y
145+
elseif ndims(Y) == 1
146+
return [ Y[1 + ((i >> 1) $ i)] for i = 0:length(Y)-1 ]
147+
else
148+
sz = [size(Y)...]
149+
tmp = Array(T, max(sz[region])) # storage for out-of-place permutations
150+
for d in region
151+
# somewhat ugly loops to do 1d permutations along dimension d
152+
na = prod(sz[d+1:end])
153+
n = sz[d]
154+
nb = prod(sz[1:d-1])
155+
sa = nb * n
156+
for ia = 0:na-1
157+
for ib = 1:nb
158+
i0 = ib + sa * ia
159+
for i = 0:n-1
160+
tmp[i+1] = Y[i0 + nb * ((i >> 1) $ i)]
161+
end
162+
for i = 0:n-1
163+
Y[i0 + nb * i] = tmp[i+1]
164+
end
165+
end
166+
end
167+
end
168+
return Y
169+
end
170+
end
171+
172+
# handle 1d case of strided arrays (loops in multidim case are too annoying)
173+
function ifwht{T<:fftwNumber}(X::StridedVector{T}, region)
174+
Y = ifwht_dyadic(X, region)
175+
176+
# Perform Gray-code permutation of Y (TODO: in-place?)
177+
if isempty(region)
178+
return Y
179+
else
180+
return [ Y[1 + ((i >> 1) $ i)] for i = 0:length(Y)-1 ]
181+
end
182+
end
183+
184+
############################################################################
185+
# Forward transforms (normalized by 1/N as in Matlab) and transforms
186+
# without the region argument:
187+
188+
for f in (:ifwht_natural, :ifwht_dyadic, :ifwht)
189+
g = symbol(string(f)[2:end])
190+
@eval begin
191+
$f(X) = $f(X, 1:ndims(X))
192+
$g(X) = scale!($f(X), normalization(X))
193+
$g(X,r) = scale!($f(X,r), normalization(X,r))
194+
end
195+
end
196+
197+
############################################################################
198+
# create (floating-point) Hadamard matrix, similar to Matlab function,
199+
# in natural order.
200+
201+
function hadamard(n::Integer)
202+
if power_of_two(n)
203+
return ifwht_natural(eye(n), 1)
204+
else
205+
# Punt. Matlab supports power-of-two multiples of 12 or 20,
206+
# but lots of other sizes are known as well (see e.g. Sloane's
207+
# web page http://neilsloane.com/hadamard/) and there is no
208+
# particularly efficient way to construct them other than
209+
# a look-up-table of known factors. [Given Hadamard matrices Hm
210+
# and Hn of sizes m and n, a Hadamard matrix of size mn is kron(Hm,Hn).]
211+
# It is not really clear what factors are important to support,
212+
# nor which of the (generally nonunique) Hadamard matrices to pick.
213+
throw(ArgumentError("non power-of-two Hadamard matrices aren't supported"))
214+
end
215+
end
216+
217+
############################################################################
218+
219+
end # module

‎test/hadamard.jl

+85
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
using Hadamard
2+
3+
H8 = [ 1 1 1 1 1 1 1 1
4+
1 1 1 1 -1 -1 -1 -1
5+
1 1 -1 -1 -1 -1 1 1
6+
1 1 -1 -1 1 1 -1 -1
7+
1 -1 -1 1 1 -1 -1 1
8+
1 -1 -1 1 -1 1 1 -1
9+
1 -1 1 -1 -1 1 -1 1
10+
1 -1 1 -1 1 -1 1 -1 ]
11+
12+
H8n = [ 1 1 1 1 1 1 1 1
13+
1 -1 1 -1 1 -1 1 -1
14+
1 1 -1 -1 1 1 -1 -1
15+
1 -1 -1 1 1 -1 -1 1
16+
1 1 1 1 -1 -1 -1 -1
17+
1 -1 1 -1 -1 1 -1 1
18+
1 1 -1 -1 -1 -1 1 1
19+
1 -1 -1 1 -1 1 1 -1 ]
20+
21+
H8d = [ 1 1 1 1 1 1 1 1
22+
1 1 1 1 -1 -1 -1 -1
23+
1 1 -1 -1 1 1 -1 -1
24+
1 1 -1 -1 -1 -1 1 1
25+
1 -1 1 -1 1 -1 1 -1
26+
1 -1 1 -1 -1 1 -1 1
27+
1 -1 -1 1 1 -1 -1 1
28+
1 -1 -1 1 -1 1 1 -1 ]
29+
30+
@assert ifwht(eye(8),1) == H8
31+
@assert ifwht(eye(8),2)' == H8
32+
@assert ifwht_natural(eye(8),1) == H8n
33+
@assert ifwht_natural(eye(8),2)' == H8n
34+
@assert ifwht_dyadic(eye(8),1) == H8d
35+
@assert ifwht_dyadic(eye(8),2)' == H8d
36+
37+
H32 = [ 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
38+
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
39+
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
40+
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
41+
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
42+
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
43+
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
44+
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
45+
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
46+
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
47+
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
48+
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
49+
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
50+
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
51+
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
52+
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
53+
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
54+
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
55+
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
56+
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
57+
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
58+
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
59+
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
60+
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
61+
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
62+
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
63+
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
64+
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
65+
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
66+
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
67+
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
68+
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 ]
69+
70+
@assert ifwht(eye(32),1) == H32
71+
@assert ifwht(eye(32),2)' == H32
72+
73+
X = reshape(sin([1:1024*32]), 1024,32);
74+
norminf(A) = max(abs(A))
75+
76+
for f in (:fwht, :fwht_natural, :fwht_dyadic)
77+
fi = symbol(string("i", f))
78+
@eval begin
79+
@assert norminf(X - $fi($f(X))) < 1e-14
80+
@assert norminf(X - $fi($f(X,1),1)) < 1e-14
81+
@assert norminf(X - $fi($f(X,2),2)) < 1e-14
82+
@assert norminf($f($f(X,1),2) - $f(X)) < 1e-14
83+
end
84+
end
85+

0 commit comments

Comments
 (0)