octave_packages/m/special-matrix/hadamard.m

   1 ## Copyright (C) 1993-2012 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 3 of the License, or (at
   8 ## your option) 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, see
  17 ## <http://www.gnu.org/licenses/>.
  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 (@nospell{Hn}) of size @var{n}-by-@var{n}.  The
  25 ## size @var{n} must be of the form @math{2^k * p} in which
  26 ## p is one of 1, 12, 20 or 28.  The returned matrix is normalized,
  27 ## meaning @w{@code{Hn(:,1) == 1}} and @w{@code{Hn(1,:) == 1}}.
  28 ##
  29 ## Some of the properties of Hadamard matrices are:
  30 ##
  31 ## @itemize @bullet
  32 ## @item
  33 ## @code{kron (Hm, Hn)} is a Hadamard matrix of size @var{m}-by-@var{n}.
  34 ##
  35 ## @item
  36 ## @code{Hn * Hn' = @var{n} * eye (@var{n})}.
  37 ##
  38 ## @item
  39 ## The rows of @nospell{Hn} are orthogonal.
  40 ##
  41 ## @item
  42 ## @code{det (@var{A}) <= abs (det (Hn))} for all @var{A} with
  43 ## @w{@code{abs (@var{A}(i, j)) <= 1}}.
  44 ##
  45 ## @item
  46 ## Multiplying any row or column by -1 and the matrix will remain a Hadamard
  47 ## matrix.
  48 ## @end itemize
  49 ## @seealso{compan, hankel, toeplitz}
  50 ## @end deftypefn
  51
  52
  53 ## Reference [1] contains a list of Hadamard matrices up to n=256.
  54 ## See code for h28 in hadamard.m for an example of how to extend
  55 ## this function for additional p.
  56 ##
  57 ## References:
  58 ## [1] A Library of Hadamard Matrices, N. J. A. Sloane
  59 ##     http://www.research.att.com/~njas/hadamard/
  60
  61 function h = hadamard (n)
  62
  63   if (nargin != 1)
  64     print_usage ();
  65   endif
  66
  67   ## Find k if n = 2^k*p.
  68   k = 0;
  69   while (n > 1 && fix (n/2) == n/2)
  70     k++;
  71     n /= 2;
  72   endwhile
  73
  74   ## Find base hadamard.
  75   ## Except for n=2^k, need a multiple of 4.
  76   if (n != 1)
  77     k -= 2;
  78   endif
  79
  80   ## Trigger error if not a multiple of 4.
  81   if (k < 0)
  82     n =- 1;
  83   endif
  84
  85   switch (n)
  86     case 1
  87       h = 1;
  88     case 3
  89       h = h12 ();
  90     case 5
  91       h = h20 ();
  92     case 7
  93       h = h28 ();
  94     otherwise
  95       error ("hadamard: N must be 2^k*p, for p = 1, 12, 20 or 28");
  96   endswitch
  97
  98   ## Build H(2^k*n) from kron(H(2^k),H(n)).
  99   h2 = [1,1;1,-1];
 100   while (true)
 101     if (fix (k/2) != k/2)
 102       h = kron (h2, h);
 103     endif
 104     k = fix (k/2);
 105     if (k == 0)
 106       break;
 107     endif
 108     h2 = kron (h2, h2);
 109   endwhile
 110
 111 endfunction
 112
 113 function h = h12 ()
 114   tu = [-1,+1,-1,+1,+1,+1,-1,-1,-1,+1,-1];
 115   tl = [-1,-1,+1,-1,-1,-1,+1,+1,+1,-1,+1];
 116   ## Note: assert (tu(2:end), tl(end:-1:2)).
 117   h = ones (12);
 118   h(2:end,2:end) = toeplitz (tu, tl);
 119 endfunction
 120
 121 function h = h20 ()
 122   tu = [+1,-1,-1,+1,+1,+1,+1,-1,+1,-1,+1,-1,-1,-1,-1,+1,+1,-1,-1];
 123   tl = [+1,-1,-1,+1,+1,-1,-1,-1,-1,+1,-1,+1,-1,+1,+1,+1,+1,-1,-1];
 124   ## Note: assert (tu(2:end), tl(end:-1:2)).
 125   h = ones (20);
 126   h(2:end,2:end) = fliplr (toeplitz (tu, tl));
 127 endfunction
 128
 129 function h = h28 ()
 130   ## Williamson matrix construction from
 131   ## http://www.research.att.com/~njas/hadamard/had.28.will.txt
 132   ## Normalized so that each row and column starts with +1
 133   h = [1 1  1  1  1  1  1  1  1 1  1  1  1 1 1 1 1 1  1 1 1 1 1  1 1  1 1  1
 134        1 1 -1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 1 1 1 1 1 -1 1 1 1 1 -1 1 -1 1 -1
 135        1 -1 1 -1 -1 -1 -1 1 -1 1 1 -1 -1 1 -1 -1 -1 -1 1 1 -1 1 -1 1 1 1 1 1
 136        1 -1 -1 1 -1 -1 -1 1 1 1 1 1 -1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 1 -1 1
 137        1 -1 -1 -1 1 -1 -1 1 1 -1 1 1 1 1 1 -1 -1 -1 1 1 1 -1 1 -1 1 -1 -1 -1
 138        1 -1 -1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1 -1 1 1 1 1 1 1 -1
 139        1 -1 -1 -1 -1 -1 1 -1 1 -1 -1 -1 1 -1 1 1 1 -1 1 1 1 1 -1 1 -1 1 -1 1
 140        1 -1 1 1 1 1 -1 -1 1 -1 -1 -1 -1 1 1 1 -1 -1 -1 -1 1 -1 -1 1 1 1 1 -1
 141        1 -1 -1 1 1 1 1 1 -1 -1 -1 -1 -1 1 -1 -1 -1 1 -1 1 1 1 1 -1 -1 1 -1 1
 142        1 1 1 1 -1 -1 -1 -1 -1 -1 1 1 1 -1 1 -1 -1 -1 -1 -1 1 1 1 -1 -1 1 1 1
 143        1 -1 1 1 1 -1 -1 -1 -1 1 -1 1 1 -1 -1 -1 1 1 -1 1 1 -1 -1 1 1 -1 -1 1
 144        1 -1 -1 1 1 1 -1 -1 -1 1 1 -1 1 -1 -1 1 1 -1 1 1 -1 -1 1 -1 -1 1 1 -1
 145        1 -1 -1 -1 1 1 1 -1 -1 1 1 1 -1 -1 1 1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 1
 146        1 1 1 1 1 -1 -1 1 1 -1 -1 -1 -1 -1 -1 1 1 -1 1 -1 -1 1 1 -1 1 -1 -1 1
 147        1 1 -1 1 1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 -1 1 -1 1 -1 1 -1 -1 1 -1
 148        1 1 -1 1 -1 -1 1 -1 1 1 1 -1 -1 -1 -1 -1 -1 1 1 -1 1 -1 1 1 1 1 -1 -1
 149        1 1 -1 1 -1 1 1 1 1 1 -1 -1 1 -1 1 -1 -1 -1 -1 1 -1 -1 -1 -1 1 -1 1 1
 150        1 1 -1 1 -1 1 -1 1 -1 1 -1 1 1 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 -1 -1
 151        1 -1 1 -1 1 -1 1 1 1 1 1 -1 1 -1 -1 1 -1 1 -1 -1 1 1 -1 -1 -1 -1 1 -1
 152        1 1 1 -1 1 -1 1 1 -1 1 -1 -1 1 1 1 -1 1 -1 -1 -1 -1 -1 1 1 -1 1 -1 -1
 153        1 1 -1 -1 1 -1 1 -1 -1 -1 -1 1 1 1 -1 1 -1 1 1 -1 -1 -1 -1 -1 1 1 1 1
 154        1 1 1 -1 -1 1 1 1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1
 155        1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 -1 -1 1 -1 1 1 1 -1 1 -1 -1 1 -1 -1 1 1
 156        1 -1 1 -1 -1 1 1 -1 1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 -1 -1 -1 1 1
 157        1 1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 -1 1 1 1 -1 1 -1 -1 -1 1 -1 -1
 158        1 -1 1 1 -1 1 1 -1 -1 -1 1 -1 1 1 1 -1 1 1 1 -1 -1 1 -1 -1 1 -1 -1 -1
 159        1 1 1 -1 -1 1 -1 -1 1 -1 1 -1 1 1 -1 1 -1 1 -1 1 -1 -1 1 1 -1 -1 -1 1
 160        1 -1 1 1 -1 -1 1 1 -1 -1 -1 1 -1 -1 1 1 -1 1 1 1 -1 -1 1 1 -1 -1 1 -1];
 161 endfunction
 162
 163
 164 %!assert (hadamard (1), 1)
 165 %!assert (hadamard (2), [1,1;1,-1])
 166 %!test
 167 %! for n = [1,2,4,8,12,24,48,20,28,2^9]
 168 %!   h = hadamard(n);
 169 %!   assert (norm (h*h' - n*eye (n)), 0);
 170 %! endfor
 171
 172 %!error hadamard ()
 173 %!error hadamard (1,2)
 174 %!error <N must be 2\^k\*p> hadamard (5)
 175