シルベスターの方法で生成したアダマール行列からウォルシュ型アダマール行列を得る方法 [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という前提．