X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fm%2Fspecial-matrix%2Fhadamard.m;fp=octave_packages%2Fm%2Fspecial-matrix%2Fhadamard.m;h=f7426c72ab258fdfad4ecfc3cac841619d270134;hp=0000000000000000000000000000000000000000;hb=1c0469ada9531828709108a4882a751d2816994a;hpb=63de9f36673d49121015e3695f2c336ea92bc278
diff --git a/octave_packages/m/special-matrix/hadamard.m b/octave_packages/m/special-matrix/hadamard.m
new file mode 100644
index 0000000..f7426c7
--- /dev/null
+++ b/octave_packages/m/special-matrix/hadamard.m
@@ -0,0 +1,175 @@
+## Copyright (C) 1993-2012 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING. If not, see
+## .
+##
+## Original version by Paul Kienzle distributed as free software in the
+## public domain.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} hadamard (@var{n})
+## Construct a Hadamard matrix (@nospell{Hn}) of size @var{n}-by-@var{n}. The
+## size @var{n} must be of the form @math{2^k * p} in which
+## p is one of 1, 12, 20 or 28. The returned matrix is normalized,
+## meaning @w{@code{Hn(:,1) == 1}} and @w{@code{Hn(1,:) == 1}}.
+##
+## Some of the properties of Hadamard matrices are:
+##
+## @itemize @bullet
+## @item
+## @code{kron (Hm, Hn)} is a Hadamard matrix of size @var{m}-by-@var{n}.
+##
+## @item
+## @code{Hn * Hn' = @var{n} * eye (@var{n})}.
+##
+## @item
+## The rows of @nospell{Hn} are orthogonal.
+##
+## @item
+## @code{det (@var{A}) <= abs (det (Hn))} for all @var{A} with
+## @w{@code{abs (@var{A}(i, j)) <= 1}}.
+##
+## @item
+## Multiplying any row or column by -1 and the matrix will remain a Hadamard
+## matrix.
+## @end itemize
+## @seealso{compan, hankel, toeplitz}
+## @end deftypefn
+
+
+## Reference [1] contains a list of Hadamard matrices up to n=256.
+## See code for h28 in hadamard.m for an example of how to extend
+## this function for additional p.
+##
+## References:
+## [1] A Library of Hadamard Matrices, N. J. A. Sloane
+## http://www.research.att.com/~njas/hadamard/
+
+function h = hadamard (n)
+
+ if (nargin != 1)
+ print_usage ();
+ endif
+
+ ## Find k if n = 2^k*p.
+ k = 0;
+ while (n > 1 && fix (n/2) == n/2)
+ k++;
+ n /= 2;
+ endwhile
+
+ ## Find base hadamard.
+ ## Except for n=2^k, need a multiple of 4.
+ if (n != 1)
+ k -= 2;
+ endif
+
+ ## Trigger error if not a multiple of 4.
+ if (k < 0)
+ n =- 1;
+ endif
+
+ switch (n)
+ case 1
+ h = 1;
+ case 3
+ h = h12 ();
+ case 5
+ h = h20 ();
+ case 7
+ h = h28 ();
+ otherwise
+ error ("hadamard: N must be 2^k*p, for p = 1, 12, 20 or 28");
+ endswitch
+
+ ## Build H(2^k*n) from kron(H(2^k),H(n)).
+ h2 = [1,1;1,-1];
+ while (true)
+ if (fix (k/2) != k/2)
+ h = kron (h2, h);
+ endif
+ k = fix (k/2);
+ if (k == 0)
+ break;
+ endif
+ h2 = kron (h2, h2);
+ endwhile
+
+endfunction
+
+function h = h12 ()
+ tu = [-1,+1,-1,+1,+1,+1,-1,-1,-1,+1,-1];
+ tl = [-1,-1,+1,-1,-1,-1,+1,+1,+1,-1,+1];
+ ## Note: assert (tu(2:end), tl(end:-1:2)).
+ h = ones (12);
+ h(2:end,2:end) = toeplitz (tu, tl);
+endfunction
+
+function h = h20 ()
+ tu = [+1,-1,-1,+1,+1,+1,+1,-1,+1,-1,+1,-1,-1,-1,-1,+1,+1,-1,-1];
+ tl = [+1,-1,-1,+1,+1,-1,-1,-1,-1,+1,-1,+1,-1,+1,+1,+1,+1,-1,-1];
+ ## Note: assert (tu(2:end), tl(end:-1:2)).
+ h = ones (20);
+ h(2:end,2:end) = fliplr (toeplitz (tu, tl));
+endfunction
+
+function h = h28 ()
+ ## Williamson matrix construction from
+ ## http://www.research.att.com/~njas/hadamard/had.28.will.txt
+ ## Normalized so that each row and column starts with +1
+ 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
+ 1 1 -1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 1 1 1 1 1 -1 1 1 1 1 -1 1 -1 1 -1
+ 1 -1 1 -1 -1 -1 -1 1 -1 1 1 -1 -1 1 -1 -1 -1 -1 1 1 -1 1 -1 1 1 1 1 1
+ 1 -1 -1 1 -1 -1 -1 1 1 1 1 1 -1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 1 -1 1
+ 1 -1 -1 -1 1 -1 -1 1 1 -1 1 1 1 1 1 -1 -1 -1 1 1 1 -1 1 -1 1 -1 -1 -1
+ 1 -1 -1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1 -1 1 1 1 1 1 1 -1
+ 1 -1 -1 -1 -1 -1 1 -1 1 -1 -1 -1 1 -1 1 1 1 -1 1 1 1 1 -1 1 -1 1 -1 1
+ 1 -1 1 1 1 1 -1 -1 1 -1 -1 -1 -1 1 1 1 -1 -1 -1 -1 1 -1 -1 1 1 1 1 -1
+ 1 -1 -1 1 1 1 1 1 -1 -1 -1 -1 -1 1 -1 -1 -1 1 -1 1 1 1 1 -1 -1 1 -1 1
+ 1 1 1 1 -1 -1 -1 -1 -1 -1 1 1 1 -1 1 -1 -1 -1 -1 -1 1 1 1 -1 -1 1 1 1
+ 1 -1 1 1 1 -1 -1 -1 -1 1 -1 1 1 -1 -1 -1 1 1 -1 1 1 -1 -1 1 1 -1 -1 1
+ 1 -1 -1 1 1 1 -1 -1 -1 1 1 -1 1 -1 -1 1 1 -1 1 1 -1 -1 1 -1 -1 1 1 -1
+ 1 -1 -1 -1 1 1 1 -1 -1 1 1 1 -1 -1 1 1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 1
+ 1 1 1 1 1 -1 -1 1 1 -1 -1 -1 -1 -1 -1 1 1 -1 1 -1 -1 1 1 -1 1 -1 -1 1
+ 1 1 -1 1 1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 -1 1 -1 1 -1 1 -1 -1 1 -1
+ 1 1 -1 1 -1 -1 1 -1 1 1 1 -1 -1 -1 -1 -1 -1 1 1 -1 1 -1 1 1 1 1 -1 -1
+ 1 1 -1 1 -1 1 1 1 1 1 -1 -1 1 -1 1 -1 -1 -1 -1 1 -1 -1 -1 -1 1 -1 1 1
+ 1 1 -1 1 -1 1 -1 1 -1 1 -1 1 1 1 -1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 -1 -1
+ 1 -1 1 -1 1 -1 1 1 1 1 1 -1 1 -1 -1 1 -1 1 -1 -1 1 1 -1 -1 -1 -1 1 -1
+ 1 1 1 -1 1 -1 1 1 -1 1 -1 -1 1 1 1 -1 1 -1 -1 -1 -1 -1 1 1 -1 1 -1 -1
+ 1 1 -1 -1 1 -1 1 -1 -1 -1 -1 1 1 1 -1 1 -1 1 1 -1 -1 -1 -1 -1 1 1 1 1
+ 1 1 1 -1 -1 1 1 1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1
+ 1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 -1 -1 1 -1 1 1 1 -1 1 -1 -1 1 -1 -1 1 1
+ 1 -1 1 -1 -1 1 1 -1 1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 -1 -1 -1 1 1
+ 1 1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 -1 1 1 1 -1 1 -1 -1 -1 1 -1 -1
+ 1 -1 1 1 -1 1 1 -1 -1 -1 1 -1 1 1 1 -1 1 1 1 -1 -1 1 -1 -1 1 -1 -1 -1
+ 1 1 1 -1 -1 1 -1 -1 1 -1 1 -1 1 1 -1 1 -1 1 -1 1 -1 -1 1 1 -1 -1 -1 1
+ 1 -1 1 1 -1 -1 1 1 -1 -1 -1 1 -1 -1 1 1 -1 1 1 1 -1 -1 1 1 -1 -1 1 -1];
+endfunction
+
+
+%!assert (hadamard (1), 1)
+%!assert (hadamard (2), [1,1;1,-1])
+%!test
+%! for n = [1,2,4,8,12,24,48,20,28,2^9]
+%! h = hadamard(n);
+%! assert (norm (h*h' - n*eye (n)), 0);
+%! endfor
+
+%!error hadamard ()
+%!error hadamard (1,2)
+%!error hadamard (5)
+