|
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 |
|
|
7 ## the Free Software Foundation; either version 2, or (at your option) |
|
|
8 ## any later version. |
|
|
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 |
|
|
16 ## along with Octave; see the file COPYING. If not, write to the Free |
|
|
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
18 ## 02110-1301, USA. |
|
|
19 ## |
|
|
20 ## Original version by Paul Kienzle distributed as free software in the |
|
|
21 ## public domain. |
|
|
22 |
|
|
23 ## -*- texinfo -*- |
|
|
24 ## @deftypefn {Function File} {} hadamard (@var{n}) |
|
|
25 ## |
|
|
26 ## Construct a Hadamard matrix @var{Hn} of size @var{n}-by-@var{n}. The |
|
|
27 ## size @var{n} must be of the form @code{2 ^ @var{k} * @var{p}} where |
|
|
28 ## @var{p} is one of 1, 12, 20 or 28. The returned matrix is normalized, |
|
|
29 ## meaning @var{Hn (:, 1) == 1} and @var{H (1, :) == 1}. |
|
|
30 ## |
|
|
31 ## Some of the properties of Hadamard matrices are: |
|
|
32 ## |
|
|
33 ## @itemize @bullet |
|
|
34 ## @item |
|
|
35 ## @code{kron (@var{Hm}, @var{Hn})} is a Hadamard matrix of size |
|
|
36 ## @var{m}-by-@var{n}. |
|
|
37 ## @item |
|
|
38 ## @code{Hn * Hn' == @var{n) * eye (@var{n})}. |
|
|
39 ## @item |
|
|
40 ## The rows of @var{Hn} are orthogonal. |
|
|
41 ## @item |
|
|
42 ## @code{det (@var{A}) <= det (@var{Hn})} for all @var{A} with |
|
|
43 ## @code{abs (@var{A} (@var{i}, @var{j})) <= 1}. |
|
|
44 ## @item |
|
|
45 ## Multiply any row or column by -1 and still have a Hadamard matrix. |
|
|
46 ## @end itemize |
|
|
47 ## |
|
|
48 ## @end deftypefn |
|
|
49 |
|
|
50 |
|
|
51 ## Reference [1] contains a list of Hadamard matrices up to n=256. |
|
|
52 ## See code for h28 in hadamard.m for an example of how to extend |
|
|
53 ## this function for additional p. |
|
|
54 ## |
|
|
55 ## References: |
|
|
56 ## [1] A Library of Hadamard Matrices, N. J. A. Sloane |
|
|
57 ## http://www.research.att.com/~njas/hadamard/ |
|
|
58 |
|
|
59 function h = hadamard (n) |
|
|
60 |
|
|
61 if (nargin != 1) |
|
|
62 print_usage (); |
|
|
63 endif |
|
|
64 |
|
|
65 ## find k if n = 2^k*p |
|
|
66 k = 0; |
|
|
67 while (n > 1 && floor (n/2) == n/2) |
|
|
68 k++; |
|
|
69 n = n/2; |
|
|
70 endwhile |
|
|
71 |
|
|
72 ## find base hadamard |
|
|
73 ## except for n=2^k, need a multiple of 4 |
|
|
74 if (n != 1) |
|
|
75 k -= 2; |
|
|
76 endif |
|
|
77 |
|
|
78 ## trigger error if not a multiple of 4 |
|
|
79 if (k < 0) |
|
|
80 n =- 1; |
|
|
81 endif |
|
|
82 |
|
|
83 switch (n) |
|
|
84 case 1 |
|
|
85 h = 1; |
|
|
86 case 3 |
|
|
87 h = h12 (); |
|
|
88 case 5 |
|
|
89 h = h20 (); |
|
|
90 case 7 |
|
|
91 h = hnormalize (h28 ()); |
|
|
92 otherwise |
|
|
93 error ("n must be 2^k*p, for p = 1, 12, 20 or 28"); |
|
|
94 endswitch |
|
|
95 |
|
|
96 ## build H(2^k*n) from kron(H(2^k),H(n)) |
|
|
97 h2 = [1,1;1,-1]; |
|
|
98 while (true) |
|
|
99 if (floor (k/2) != k/2) |
|
|
100 h = kron (h2, h); |
|
|
101 endif |
|
|
102 k = floor (k/2); |
|
|
103 if (k == 0) |
|
|
104 break; |
|
|
105 endif |
|
|
106 h2 = kron (h2, h2); |
|
|
107 endwhile |
|
|
108 endfunction |
|
|
109 |
|
|
110 function h = h12 () |
|
|
111 tu=[-1,+1,-1,+1,+1,+1,-1,-1,-1,+1,-1]; |
|
|
112 tl=[-1,-1,+1,-1,-1,-1,+1,+1,+1,-1,+1]; |
|
|
113 ## note: assert(tu(2:end),tl(end:-1:2)) |
|
|
114 h = ones (12); |
|
|
115 h(2:end,2:end) = toeplitz (tu, tl); |
|
|
116 endfunction |
|
|
117 |
|
|
118 |
|
|
119 function h = h20 () |
|
|
120 tu=[+1,-1,-1,+1,+1,+1,+1,-1,+1,-1,+1,-1,-1,-1,-1,+1,+1,-1,-1]; |
|
|
121 tl=[+1,-1,-1,+1,+1,-1,-1,-1,-1,+1,-1,+1,-1,+1,+1,+1,+1,-1,-1]; |
|
|
122 ## note: assert(tu(2:end),tl(end:-1:2)) |
|
|
123 h = ones (20); |
|
|
124 h(2:end,2:end) = fliplr (toeplitz (tu, tl)); |
|
|
125 endfunction |
|
|
126 |
|
|
127 function h = hnormalize (h) |
|
|
128 ## Make sure each row and column starts with +1 |
|
|
129 h(h(:,1)==-1,:) *= -1; |
|
|
130 h(:,h(1,:)==-1) *= -1; |
|
|
131 endfunction |
|
|
132 |
|
|
133 function h = h28 () |
|
|
134 ## Williamson matrix construction from |
|
|
135 ## http://www.research.att.com/~njas/hadamard/had.28.will.txt |
|
|
136 s = ['+------++----++-+--+-+--++--'; |
|
|
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 '--++-++++-+-+----+++-----+--'; |
|
|
162 '+--++-+-++-+-+----+++-----+-'; |
|
|
163 '++--++-+-++-+-+----++------+']; |
|
5901
|
164 ## Without this, the assignment of -1 will not work properly |
|
|
165 ## (compatibility forces LHS(idx) = ANY_VAL to keep the LHS logical |
|
|
166 ## instead of widening to a type that can represent ANY_VAL). |
|
|
167 h = double (s=='+'); |
|
5827
|
168 h(!h) = -1; |
|
|
169 endfunction |
|
|
170 |
|
|
171 %!assert(hadamard(1),1) |
|
|
172 %!assert(hadamard(2),[1,1;1,-1]) |
|
|
173 %!test |
|
|
174 %! for n=[1,2,4,8,12,24,48,20,28,2^9] |
|
|
175 %! h=hadamard(n); assert(norm(h*h'-n*eye(n)),0); |
|
|
176 %! end |
