シルベスターの方法で生成したアダマール行列からウォルシュ型アダマール行列を得る方法 (C++)

信号処理や統計学で使われるアダマール行列は,以下のようにシルベスターの方法で再帰的に生成するのが(おそらく)一般的.


  \begin{align}
    \mathbf{H}_{2} &=\left[
    \begin{array}{rr}
      1 & 1 \\\
      1 & -1
    \end{array}\right]\\\
    \mathbf{H}_{2^{k}} &= \mathbf{H}_{2}\otimes\mathbf{H}_{2^{k-1}}
  \end{align}

ここで \otimesクロネッカー積を表す.例として \mathbf{H}_{8} を示す.


  \begin{equation}
  \mathbf{H}_{8} =\left[
  \begin{array}{rrrrrrrr}
    1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\\
    1 & -1 & 1 & -1 & 1 & -1 & 1 & -1\\\
    1 & 1 & -1 & -1 & 1 & 1 & -1 & -1\\\
    1 & -1 & -1 & 1 & 1 & -1 & -1 & 1\\\
    1 & 1 & 1 & 1 & -1 & -1 & -1 & -1\\\
    1 & -1 & 1 & -1 & -1 & 1 & -1 & 1\\\
    1 & 1 & -1 & -1 & -1 & -1 & 1 & 1\\\
    1 & -1 & -1 & 1 & -1 & 1 & 1 & -1
  \end{array}\right]
  \end{equation}

さて,この方法で生成したアダマール行列について,各行の符号反転の回数(周波数)を小さい順に並べ替えた方が都合がいい場合がある.そのような行列をウォルシュ型アダマール行列と呼んだりするが,次のような方法で並べ替えることができる.

/**
 * @brief ビット逆転(32 bit 用)
 * @param[in] a 32 bit 符号なし整数
 * @return ビット逆転した 32 bit 符号なし整数
 */
unsigned int bitreverse( const unsigned int a ) {
    unsigned int x = a;
    x = ( ( ( x & 0xaaaaaaaa ) >> 1 ) | ( ( x & 0x55555555 ) << 1 ) );
    x = ( ( ( x & 0xcccccccc ) >> 2 ) | ( ( x & 0x33333333 ) << 2 ) );
    x = ( ( ( x & 0xf0f0f0f0 ) >> 4 ) | ( ( x & 0x0f0f0f0f ) << 4 ) );
    x = ( ( ( x & 0xff00ff00 ) >> 8 ) | ( ( x & 0x00ff00ff ) << 8 ) );
    return( ( x >> 16 ) | ( x << 16 ) );
}

/**
 * @brief 並び替えのための配列生成
 * @param[in] k 行列の次数 (2^k)
 * @param[out] v 並び替えのための配列
 */
void permutation( const unsigned int k,
                  std::vector< unsigned int > &v ) {
    const unsigned int N = 1 << k;
    v.resize(N);
    for ( unsigned int i = 0; i < N; i++ ) {
        unsigned int j = bitreverse(i) >> ( 31 - k );
        j = j ^ ( j >> 1 );
        j = ( j << ( 32 - k ) ) >> ( 32 - k );
        v[i] = j;
    }
}

たとえばpermutation( 3, v )とした場合はv = [0 4 6 2 3 7 5 1]となるので,この順に行を並べ替えれば8次のウォルシュ型アダマール行列を得られる.上のプログラムは当然ながらsizeof( unsigned int ) == 4という前提.