Fast Sequency-ordered Walsh-Hadamard Transform

Question

EDIT You can checkout my implementation on Github: https://github.com/Sheljohn/WalshHadamard

---

I am looking for an implementation, or indications on how to implement, the sequency-ordered Fast Walsh Hadamard transform (see this and this).

I slightly adapted a very nice implementation found online:

// (a,b) -> (a+b,a-b) without overflow
void rotate( long& a, long& b )
{
    static long t;
    t = a;
    a = a + b;
    b = t - b;
}

// Integer log2
long ilog2( long x )
{
    long l2 = 0;
    for (; x; x >>=1) ++l2;
    return l2;
}

/**
 * Fast Walsh-Hadamard transform
 */
void fwht( std::vector<long>& data )
{
    const long l2 = ilog2(data.size()) - 1;

    for (long i = 0; i < l2; ++i)
    {
        for (long j = 0; j < (1 << l2); j += 1 << (i+1))
        for (long k = 0; k < (1 << i ); ++k)
            rotate( data[j + k], data[j + k + (1<<i)] );
    }
}

but it does not compute the WHT in sequency order (the natural Hadamard matrix is used implicitly). Note that in the code above (and if you try it), the size of data needs to be a power of 2.

My question is: is there a simple adaptation of this implementation that gives the sequency-ordered FWHT?

A possible solution would be to write a small function to compute dynamically the elements of Hn (the Hadamard matrix of order n), count the number of zero crossings, and create a ranking of the rows, but I am wondering whether there is a smarter way. Thanks in advance for any input! Cheers

Answer 1

As indicated here (linked from within your reference):

The sequency ordering of the rows of the Walsh matrix can be derived from the ordering of the Hadamard matrix by first applying the bit-reversal permutation and then the Gray code permutation.

There are various implementations of bit-reversal algorithm such as this:

// Bit-reversal
// adapted from http://www.idi.ntnu.no/~elster/pubs/elster-bit-rev-1989.pdf
void bitrev(int t, std::vector<long>& c)
{
  long n = 1<<t;
  long L = 1;
  c[0]  = 0;
  for (int q=0; q<t; ++q)
  {
    n /= 2;
    for (long j=0; j<L; ++j)
    {
      c[L+j] = c[j] + n;
    }
    L *= 2;
  }
}

The gray code can be obtained from here:

/*
    The purpose of this function is to convert an unsigned
    binary number to reflected binary Gray code.

    The operator >> is shift right. The operator ^ is exclusive or.
*/
unsigned int binaryToGray(unsigned int num)
{
    return (num >> 1) ^ num;
}

These can be combined to yields the final permutation:

// Compute a permutation of size 2^order 
// to reorder the Fast Walsh-Hadamard transform's output 
// into the Walsh-ordered (sequency-ordered)
void sequency_permutation(long order, std::vector<long>& p)
{
  long n = 1<<order;

  std::vector<long> tmp(n);
  bitrev(order, tmp);
  p.resize(n);
  for (long i=0; i<n; ++i)
  {
    p[i] = tmp[binaryToGray(i)];
  }
}

All that's left to do is to apply the permutation to the normal Walsh-Hadamard Transform output.

void permuted_fwht(std::vector<long>& data, const std::vector<long>& permutation)
{
  std::vector<long> tmp = data;
  fwht(tmp);

  for (long i=0; i<data.size(); ++i)
  {
    data[i] = tmp[permutation[i]];
  }
}

Note that the permutation is fixed for a given data size, so it only needs to be computed once (assuming you are processing multiple blocks of data). So, putting it all together you would get something such as:

  std::vector<long> p;

  const long order = ilog2(data_block_size) - 1;
  sequency_permutation(order, p);

  permuted_fwht( data_block_1, p);
  permuted_fwht( data_block_2, p);
  //...