|
5827
|
1 ## Copyright (C) 2004 Paul Kienzle |
|
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
|
6 ## under the terms of the GNU General Public License as published by |
|
7016
|
7 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
8 ## your option) any later version. |
|
5827
|
9 ## |
|
|
10 ## Octave is distributed in the hope that it will be useful, but |
|
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
13 ## General Public License for more details. |
|
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
7016
|
16 ## along with Octave; see the file COPYING. If not, see |
|
|
17 ## <http://www.gnu.org/licenses/>. |
|
5827
|
18 ## |
|
|
19 ## Original version by Paul Kienzle distributed as free software in the |
|
|
20 ## public domain. |
|
|
21 |
|
|
22 ## -*- texinfo -*- |
|
|
23 ## @deftypefn {Function File} {} hadamard (@var{n}) |
|
|
24 ## Construct a Hadamard matrix @var{Hn} of size @var{n}-by-@var{n}. The |
|
6516
|
25 ## size @var{n} must be of the form @code{2 ^ @var{k} * @var{p}} in which |
|
5827
|
26 ## @var{p} is one of 1, 12, 20 or 28. The returned matrix is normalized, |
|
6516
|
27 ## meaning @code{Hn(:,1) == 1} and @code{H(1,:) == 1}. |
|
5827
|
28 ## |
|
|
29 ## Some of the properties of Hadamard matrices are: |
|
|
30 ## |
|
|
31 ## @itemize @bullet |
|
|
32 ## @item |
|
|
33 ## @code{kron (@var{Hm}, @var{Hn})} is a Hadamard matrix of size |
|
|
34 ## @var{m}-by-@var{n}. |
|
|
35 ## @item |
|
6248
|
36 ## @code{Hn * Hn' == @var{n} * eye (@var{n})}. |
|
5827
|
37 ## @item |
|
|
38 ## The rows of @var{Hn} are orthogonal. |
|
|
39 ## @item |
|
|
40 ## @code{det (@var{A}) <= det (@var{Hn})} for all @var{A} with |
|
|
41 ## @code{abs (@var{A} (@var{i}, @var{j})) <= 1}. |
|
|
42 ## @item |
|
|
43 ## Multiply any row or column by -1 and still have a Hadamard matrix. |
|
|
44 ## @end itemize |
|
|
45 ## |
|
|
46 ## @end deftypefn |
|
|
47 |
|
|
48 |
|
|
49 ## Reference [1] contains a list of Hadamard matrices up to n=256. |
|
|
50 ## See code for h28 in hadamard.m for an example of how to extend |
|
|
51 ## this function for additional p. |
|
|
52 ## |
|
|
53 ## References: |
|
|
54 ## [1] A Library of Hadamard Matrices, N. J. A. Sloane |
|
|
55 ## http://www.research.att.com/~njas/hadamard/ |
|
|
56 |
|
|
57 function h = hadamard (n) |
|
|
58 |
|
|
59 if (nargin != 1) |
|
|
60 print_usage (); |
|
|
61 endif |
|
|
62 |
|
|
63 ## find k if n = 2^k*p |
|
|
64 k = 0; |
|
|
65 while (n > 1 && floor (n/2) == n/2) |
|
|
66 k++; |
|
|
67 n = n/2; |
|
|
68 endwhile |
|
|
69 |
|
|
70 ## find base hadamard |
|
|
71 ## except for n=2^k, need a multiple of 4 |
|
|
72 if (n != 1) |
|
|
73 k -= 2; |
|
|
74 endif |
|
|
75 |
|
|
76 ## trigger error if not a multiple of 4 |
|
|
77 if (k < 0) |
|
|
78 n =- 1; |
|
|
79 endif |
|
|
80 |
|
|
81 switch (n) |
|
|
82 case 1 |
|
|
83 h = 1; |
|
|
84 case 3 |
|
|
85 h = h12 (); |
|
|
86 case 5 |
|
|
87 h = h20 (); |
|
|
88 case 7 |
|
|
89 h = hnormalize (h28 ()); |
|
|
90 otherwise |
|
|
91 error ("n must be 2^k*p, for p = 1, 12, 20 or 28"); |
|
|
92 endswitch |
|
|
93 |
|
|
94 ## build H(2^k*n) from kron(H(2^k),H(n)) |
|
|
95 h2 = [1,1;1,-1]; |
|
|
96 while (true) |
|
|
97 if (floor (k/2) != k/2) |
|
|
98 h = kron (h2, h); |
|
|
99 endif |
|
|
100 k = floor (k/2); |
|
|
101 if (k == 0) |
|
|
102 break; |
|
|
103 endif |
|
|
104 h2 = kron (h2, h2); |
|
|
105 endwhile |
|
|
106 endfunction |
|
|
107 |
|
|
108 function h = h12 () |
|
|
109 tu=[-1,+1,-1,+1,+1,+1,-1,-1,-1,+1,-1]; |
|
|
110 tl=[-1,-1,+1,-1,-1,-1,+1,+1,+1,-1,+1]; |
|
|
111 ## note: assert(tu(2:end),tl(end:-1:2)) |
|
|
112 h = ones (12); |
|
|
113 h(2:end,2:end) = toeplitz (tu, tl); |
|
|
114 endfunction |
|
|
115 |
|
|
116 |
|
|
117 function h = h20 () |
|
|
118 tu=[+1,-1,-1,+1,+1,+1,+1,-1,+1,-1,+1,-1,-1,-1,-1,+1,+1,-1,-1]; |
|
|
119 tl=[+1,-1,-1,+1,+1,-1,-1,-1,-1,+1,-1,+1,-1,+1,+1,+1,+1,-1,-1]; |
|
|
120 ## note: assert(tu(2:end),tl(end:-1:2)) |
|
|
121 h = ones (20); |
|
|
122 h(2:end,2:end) = fliplr (toeplitz (tu, tl)); |
|
|
123 endfunction |
|
|
124 |
|
|
125 function h = hnormalize (h) |
|
|
126 ## Make sure each row and column starts with +1 |
|
|
127 h(h(:,1)==-1,:) *= -1; |
|
|
128 h(:,h(1,:)==-1) *= -1; |
|
|
129 endfunction |
|
|
130 |
|
|
131 function h = h28 () |
|
|
132 ## Williamson matrix construction from |
|
|
133 ## http://www.research.att.com/~njas/hadamard/had.28.will.txt |
|
|
134 s = ['+------++----++-+--+-+--++--'; |
|
|
135 '-+-----+++-----+-+--+-+--++-'; |
|
|
136 '--+-----+++---+-+-+----+--++'; |
|
|
137 '---+-----+++---+-+-+-+--+--+'; |
|
|
138 '----+-----+++---+-+-+++--+--'; |
|
|
139 '-----+-----++++--+-+--++--+-'; |
|
|
140 '------++----++-+--+-+--++--+'; |
|
|
141 '--++++-+-------++--+++-+--+-'; |
|
|
142 '---++++-+-----+-++--+-+-+--+'; |
|
|
143 '+---+++--+----++-++--+-+-+--'; |
|
|
144 '++---++---+----++-++--+-+-+-'; |
|
|
145 '+++---+----+----++-++--+-+-+'; |
|
|
146 '++++--------+-+--++-++--+-+-'; |
|
|
147 '-++++--------+++--++--+--+-+'; |
|
|
148 '-+-++-++--++--+--------++++-'; |
|
|
149 '+-+-++--+--++--+--------++++'; |
|
|
150 '-+-+-++--+--++--+----+---+++'; |
|
|
151 '+-+-+-++--+--+---+---++---++'; |
|
|
152 '++-+-+-++--+------+--+++---+'; |
|
|
153 '-++-+-+-++--+------+-++++---'; |
|
|
154 '+-++-+---++--+------+-++++--'; |
|
|
155 '-++--++-+-++-+++----++------'; |
|
|
156 '+-++--++-+-++-+++-----+-----'; |
|
|
157 '++-++---+-+-++-+++-----+----'; |
|
|
158 '-++-++-+-+-+-+--+++-----+---'; |
|
|
159 '--++-++++-+-+----+++-----+--'; |
|
|
160 '+--++-+-++-+-+----+++-----+-'; |
|
|
161 '++--++-+-++-+-+----++------+']; |
|
5901
|
162 ## Without this, the assignment of -1 will not work properly |
|
|
163 ## (compatibility forces LHS(idx) = ANY_VAL to keep the LHS logical |
|
|
164 ## instead of widening to a type that can represent ANY_VAL). |
|
|
165 h = double (s=='+'); |
|
5827
|
166 h(!h) = -1; |
|
|
167 endfunction |
|
|
168 |
|
|
169 %!assert(hadamard(1),1) |
|
|
170 %!assert(hadamard(2),[1,1;1,-1]) |
|
|
171 %!test |
|
|
172 %! for n=[1,2,4,8,12,24,48,20,28,2^9] |
|
|
173 %! h=hadamard(n); assert(norm(h*h'-n*eye(n)),0); |
|
|
174 %! end |
