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