--- a/main/signal/NEWS	Sun Mar 10 22:32:12 2013 +0000
+++ b/main/signal/NEWS	Mon Mar 11 02:25:15 2013 +0000
@@ -7,6 +7,8 @@
  ** The following functions are new:
 
       findpeaks
+      fwht
+      ifwht
 
  ** The function `blackmanharris' was fixed to return a column vector.
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/signal/inst/fwht.m	Mon Mar 11 02:25:15 2013 +0000
@@ -0,0 +1,79 @@
+## Copyright (C) 2013 Mike Miller <mtmiller@ieee.org>
+##
+## This program 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.
+##
+## This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {} fwht (@var{x})
+## @deftypefnx {Function File} {} fwht (@var{x}, @var{n})
+## @deftypefnx {Function File} {} fwht (@var{x}, @var{n}, @var{order})
+## Compute the Walsh-Hadamard transform of @var{x} using the Fast
+## Walsh-Hadamard Transform (FWHT) algorithm.  If the input is a matrix,
+## the FWHT is calculated along the columns of @var{x}.
+##
+## The number of elements of @var{x} must be a power of 2; if not, the
+## input will be extended and filled with zeros.  If a second argument
+## is given, the input is truncated or extended to have length @var{n}.
+##
+## The third argument specifies the @var{order} in which the returned
+## Walsh-Hadamard transform coefficients should be arranged.  The
+## @var{order} may be any of the following strings:
+##
+## @table @asis
+## @item "sequency"
+## The coefficients are returned in sequency order.  This is the default
+## if @var{order} is not given.
+##
+## @item "hadamard"
+## The coefficients are returned in Hadamard order.
+##
+## @item "dyadic"
+## The coefficients are returned in Gray code order.
+## @end table
+##
+## @seealso{ifwht}
+## @end deftypefn
+
+## Author: Mike Miller
+
+function y = fwht (x, n, order)
+
+  if (nargin < 1 || nargin > 3)
+    print_usage ();
+  elseif (nargin == 1)
+    n = order = [];
+  elseif (nargin == 2)
+    order = [];
+  endif
+
+  [y, n] = __fwht_opts__ ("fwht", x, n, order);
+  y /= n;
+
+endfunction
+
+%!assert (isempty (fwht ([])));
+%!assert (fwht (zeros (16)), zeros (16));
+%!assert (fwht (ones (16, 1)), [1; (zeros (15, 1))]);
+%!assert (fwht (zeros (17, 1)), zeros (32, 1));
+%!assert (fwht ([1 -1 1 -1 1 -1 1 -1]), [0 0 0 0 0 0 0 1]);
+
+%!test
+%! x = randi (16, 16);
+%! assert (ifwht (fwht (x)), x);
+
+%% Test input validation
+%!error fwht ();
+%!error fwht (1, 2, 3, 4);
+%!error fwht (0, 0);
+%!error fwht (0, 5);
+%!error fwht (0, [], "invalid");

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/signal/inst/ifwht.m	Mon Mar 11 02:25:15 2013 +0000
@@ -0,0 +1,74 @@
+## Copyright (C) 2013 Mike Miller <mtmiller@ieee.org>
+##
+## This program 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.
+##
+## This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {} ifwht (@var{x})
+## @deftypefnx {Function File} {} ifwht (@var{x}, @var{n})
+## @deftypefnx {Function File} {} ifwht (@var{x}, @var{n}, @var{order})
+## Compute the inverse Walsh-Hadamard transform of @var{x} using the
+## Fast Walsh-Hadamard Transform (FWHT) algorithm.  If the input is a
+## matrix, the inverse FWHT is calculated along the columns of @var{x}.
+##
+## The number of elements of @var{x} must be a power of 2; if not, the
+## input will be extended and filled with zeros.  If a second argument
+## is given, the input is truncated or extended to have length @var{n}.
+##
+## The third argument specifies the @var{order} in which the returned
+## inverse Walsh-Hadamard transform coefficients should be arranged.
+## The @var{order} may be any of the following strings:
+##
+## @table @asis
+## @item "sequency"
+## The coefficients are returned in sequency order.  This is the default
+## if @var{order} is not given.
+##
+## @item "hadamard"
+## The coefficients are returned in Hadamard order.
+##
+## @item "dyadic"
+## The coefficients are returned in Gray code order.
+## @end table
+##
+## @seealso{fwht}
+## @end deftypefn
+
+## Author: Mike Miller
+
+function y = ifwht (x, n, order)
+
+  if (nargin < 1 || nargin > 3)
+    print_usage ();
+  elseif (nargin == 1)
+    n = order = [];
+  elseif (nargin == 2)
+    order = [];
+  endif
+
+  y = __fwht_opts__ ("ifwht", x, n, order);
+
+endfunction
+
+%!assert (isempty (ifwht ([])));
+%!assert (ifwht (zeros (16)), zeros (16));
+%!assert (ifwht ([1; (zeros (15, 1))]), ones (16, 1));
+%!assert (ifwht (zeros (17, 1)), zeros (32, 1));
+%!assert (ifwht ([0 0 0 0 0 0 0 1]), [1 -1 1 -1 1 -1 1 -1]);
+
+%% Test input validation
+%!error ifwht ();
+%!error ifwht (1, 2, 3, 4);
+%!error ifwht (0, 0);
+%!error ifwht (0, 5);
+%!error ifwht (0, [], "invalid");

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/signal/inst/private/__fwht_opts__.m	Mon Mar 11 02:25:15 2013 +0000
@@ -0,0 +1,93 @@
+## Copyright (C) 2013 Mike Miller <mtmiller@ieee.org>
+##
+## This program 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.
+##
+## This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{y}, @var{len}] =} __fwht_opts__ (@var{caller}, @var{x}, @var{n}, @var{order})
+## Undocumented internal function.
+## @end deftypefn
+
+## Author: Mike Miller
+
+function [y, len] = __fwht_opts__ (caller, x, n, order)
+
+  if (nargin < 1)
+    print_usage ();
+  elseif (nargin != 4)
+    print_usage (caller);
+  endif
+
+  [nr, nc] = size (x);
+  
+  if (isempty (n))
+    if (nr == 1)
+      n = 2^nextpow2 (nc);
+    else
+      n = 2^nextpow2 (nr);
+    endif
+  elseif (!(isscalar (n) && n == fix (n) && n > 0))
+    error ("%s: N must be a positive scalar", caller);
+  else
+    f = log2(n);
+    if (f != fix (f))
+      error ("%s: N must be a power of 2", caller);
+    endif
+  endif
+
+  if (isempty (order))
+    order = "sequency";
+  endif
+  if (!(strncmp (order, "dyadic", 6) ||
+        strncmp (order, "hadamard", 8) ||
+        strncmp (order, "sequency", 8)))
+    error ("%s: invalid order option", caller);
+  endif
+
+  if (nr == 1)
+    nc = n;
+    x = x(:);
+  else
+    nr = n;
+  endif
+
+  ## Zero-based index for normal Hadamard ordering
+  idx = 0:n-1;
+
+  ## Gray code permutation of index for alternate orderings
+  idx_bin = dec2bin (idx) - "0";
+  idx_bin_a = idx_bin(:,1:end-1);
+  idx_bin_b = idx_bin(:,2:end);
+  idx_bin(:,2:end) = mod (idx_bin_a + idx_bin_b, 2);
+  idx_bin = char (idx_bin + "0");
+
+  if (strncmp (order, "dyadic", 6))
+    idx = bin2dec (idx_bin) + 1;
+  elseif (strncmp (order, "sequency", 8))
+    idx = bin2dec (fliplr (idx_bin)) + 1;
+  else
+    idx += 1;
+  endif
+
+  len = n;
+  x = postpad (x, len);
+
+  if (len < 2)
+    y = x;
+  else
+    y = __fwht__ (x);
+  endif
+
+  y = reshape (y(idx,:), nr, nc);
+
+endfunction

--- a/main/signal/src/Makefile	Sun Mar 10 22:32:12 2013 +0000
+++ b/main/signal/src/Makefile	Mon Mar 11 02:25:15 2013 +0000
@@ -1,6 +1,6 @@
 MKOCTFILE ?= mkoctfile
 
-all: cl2bp.oct remez.oct medfilt1.oct sosfilt.oct upfirdn.oct
+all: __fwht__.oct cl2bp.oct remez.oct medfilt1.oct sosfilt.oct upfirdn.oct
 
 %.oct: %.cc
 	$(MKOCTFILE) -Wall -s $<

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/signal/src/__fwht__.cc	Mon Mar 11 02:25:15 2013 +0000
@@ -0,0 +1,62 @@
+// Copyright (C) 2013 Mike Miller <mtmiller@ieee.org>
+//
+// This program 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.
+//
+// This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#include <octave/oct.h>
+
+template <typename T>
+static void do_fwht (T *x, octave_idx_type n)
+{
+  double *xa = &x[0];
+  double *xb = &x[n/2];
+
+  // Perform the FWHT butterfly on the input
+  for (octave_idx_type i = 0; i < n/2; i++)
+    {
+      double ya = xa[i] + xb[i];
+      double yb = xa[i] - xb[i];
+      xa[i] = ya;
+      xb[i] = yb;
+    }
+
+  // Divide and conquer
+  if (n > 2)
+    {
+      do_fwht (xa, n/2);
+      do_fwht (xb, n/2);
+    }
+}
+
+DEFUN_DLD (__fwht__, args, ,
+  "-*- texinfo -*-\n"
+"@deftypefn {Loadable Function} {} __fwht__ (@var{x})\n"
+"Undocumented internal function.\n"
+"@end deftypefn\n")
+{
+  int nargin = args.length ();
+  if (nargin != 1)
+    print_usage ();
+  else
+    {
+      NDArray x = args(0).array_value ();
+      octave_idx_type n = x.rows ();
+      double *data = x.fortran_vec ();
+      for (octave_idx_type i = 0; i < x.columns (); i++)
+        {
+          do_fwht (&data[i*n], n);
+        }
+      return octave_value (x);
+    }
+  return octave_value ();
+}