function x = ffwt ( n, x )

%*****************************************************************************80
%
%% ffwt() performs an in-place fast Walsh transform.
%
%  Discussion:
%
%    This routine performs a fast Walsh transform on an input series X
%    leaving the transformed results in X.
%    X is dimensioned N, which must be a power of 2.
%    The results of this Walsh transform are in sequency order.
%
%    The output sequence could be normalized by dividing by N.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    16 March 2011
%
%  Author:
%
%    Ken Beauchamp
%
%  Reference:
%
%    Ken Beauchamp,
%    Walsh functions and their applications,
%    Academic Press, 1975,
%    ISBN: 0-12-084050-2,
%    LC: QA404.5.B33.
%
%  Input:
%
%    integer N, the number of items in X.
%    N must be a power of 2.
%
%    real X(N), the data to be transformed.
%
%  Output:
%
%    real X(N), the transformed data.
%
  m = i4_log_2 ( n );

  for i = 1 : m
    two_power(i) = 2^( m - i );
  end

  for l = 1 : m

    nz = 2^( l - 1 );
    nzi = 2 * nz;
    nzn = floor ( n / nzi );
    nz2 = floor ( nz / 2 );
    if ( nz2 == 0 )
      nz2 = 1;
    end

    for i = 1 : nzn

      js = ( i - 1 ) * nzi;
      z = 1.0;

      for ii = 1 : 2

        for j = 1 : nz2
          js = js + 1;
          j2 = js + nz;
          hold = x(js) + z * x(j2);
          z = - z;
          x(j2) = x(js) + z * x(j2);
          x(js) = hold;
          z = - z;
        end

        if ( l == 1 )
          break
        end

        z = - 1.0;

      end

    end
  end
%
%  Bit reversal section.
%
  nw = 0;
  for k = 1 : n
%
%  Choose correct index and switch elements if not already switched.
%
    if ( k < nw + 1 )
      hold = x(nw+1);
      x(nw+1) = x(k);
      x(k) = hold;
    end
%
%  Bump up series by 1.
%
    for i = 1 : m

      ii = i;
      if ( nw < two_power(i) )
        break
      end
      mw = floor ( nw / two_power(i) );
      mw1 = floor ( mw / 2 );
      if ( mw <= 2 * mw1 )
        break
      end

      nw = nw - two_power(i);

    end

    nw = nw + two_power(ii);

  end

  return
end
 
