MDCTで使う窓関数の逆変換

MDCTの窓関数の逆変換を実装してみた。

窓関数処理

    for (int i = 0; i < (nsmpl >> 1); i++) {
        p_work[i] = -p_x[3 * nsmpl / 2 - 1 - i] * p_w[nsmpl / 2 + i] -
                         p_x[3 * nsmpl / 2 + i] * p_w[nsmpl / 2 - 1 - i];
        p_work[i]] = p_x[i] * p_w[i] - p_x[nsmpl - 1 - i] * p_w[nsmpl - 1 - i];
    }

というMDCTの窓関数処理の逆変換を実装してみました。

• nsmpl サンプル数
• p_x 入力列の配列　サイズはnsmpl*2
• p_w 窓関数の配列 サイズはnsmpl

　サンプル数４の場合

まず、nsmpl を 4 としてループを展開してみる。

      p_work[0] = -p_x[5] * p_w[2] - p_x[6] * p_w[1];
      p_work[2] =  p_x[0] * p_w[0] - p_x[3] * p_w[3];
      p_work[1] = -p_x[4] * p_w[3] - p_x[7] * p_w[0];
      p_work[3] =  p_x[1] * p_w[1] - p_x[2] * p_w[2];

MDCTなので窓を半分ずつ重複させながら変換する。

この窓が半分ずつ重なっている部分の情報を使うと逆変換でp_x[]を取り出すことができる。

一つ前の変換のp_work[]をp_prv_work[]、p_x[]をp_prv_x[]とする

      p_prv_work[0] = -p_prv_x[5] * p_w[2] - p_prv_x[6] * p_w[1];

重なり部分の要素数はnsmpl == 4 の場合は重なりも4となる。
そうすると、一つ前のp_prv_x[5] と 現在のp_x[1] 、一つ前のp_prv_x[6]と 現在のp_x[2]が同一の値となる。

現在の情報

     sig_cur[0] = p_work[2];
     sig_cur[1] = p_work[3];
     sig_cur[2] = p_work[3];
     sig_cur[3] = p_work[p];

一つ前の情報

     sig_prv[0] = p_prv_work[1];
     sig_prv[1] = p_prv_work[0];
     sig_prv[2] = -p_prv_work[0];
     sig_prv[3] = -p_prv_work[1];

として以下のような変換式が得られる。

     sig_cur[0] * p_w[0] - sig_prv[0] * p_w[3]
　　 = p_work[2] * p_w[0] - p_prv_work[1] * p_w[3]
　　 = (p_x[0] * p_w[0] - p_x[3] * p_w[3]) * p_w[0] -  (-p_prv_x[4] * p_w[3] - p_prv_x[7] * p_w[0];)* p_w[3]
　　 = (p_x[0] * p_w[0] - p_x[3] * p_w[3]) * p_w[0] -  (-p_x[0] * p_w[3] - p_x[3] * p_w[0];)* p_w[3]
     = p_x[0] * p_w[0] * p_w[0] - p_x[3] * p_w[3] * p_w[0] + p_x[0] * p_w[3] * p_w[3] + p_x[3] * p_w[3] * p_w[0]
     = p_x[0] * ( p_w[0] * p_w[0] + p_w[3] * p_w[3]) 

( p_w[0] * p_w[0] + p_w[3] * p_w[3])は窓関数の配列から計算で導出できるので、割り算するとp_x[0]が得られる。

同様の計算でp_x[1]、p_x[2]、p_x[3]が得られる。

    p_x[0] : (sig_cur[0] * p_w[0] - sig_prv[0] * p_w[3]) / (p_w[0] * p_w[0] + p_w[3] * p_w[3]);
    p_x[1] : (sig_cur[1] * p_w[1] - sig_prv[1] * p_w[2]) / (p_w[1] * p_w[1] + p_w[2] * p_w[2]);
    p_x[2] : (sig_cur[2] * p_w[2] - sig_prv[2] * p_w[1]) / (-p_w[1] * p_w[1] - p_w[2] * p_w[2]);
    p_x[3] : (sig_cur[3] * p_w[3] - sig_prv[3] * p_w[0]) / (-p_w[0] * p_w[0] - p_w[3] * p_w[3]);

サンプル数８の場合

サンプル数８の場合の変換式は以下となる。

    sig_cur[0] = p_work[4];
    sig_cur[1] = p_work[5];
    sig_cur[2] = p_work[6];
    sig_cur[3] = p_work[7];
    sig_cur[4] = p_work[7];
    sig_cur[5] = p_work[6];
    sig_cur[6] = p_work[5];
    sig_cur[7] = p_work[4];

    p_z[0] = (sig_cur[0] * p_w[0] - sig_prv[0] * p_w[7]) / (p_w[0] * p_w[0] + p_w[7] * p_w[7]);
    p_z[1] = (sig_cur[1] * p_w[1] - sig_prv[1] * p_w[6]) / (p_w[1] * p_w[1] + p_w[6] * p_w[6]);
    p_z[2] = (sig_cur[2] * p_w[2] - sig_prv[2] * p_w[5]) / (p_w[2] * p_w[2] + p_w[5] * p_w[5]);
    p_z[3] = (sig_cur[3] * p_w[3] - sig_prv[3] * p_w[4]) / (p_w[3] * p_w[3] + p_w[4] * p_w[4]);
    p_z[4] = (sig_cur[4] * p_w[4] - sig_prv[4] * p_w[3]) / (-p_w[3] * p_w[3] - p_w[4] * p_w[4]);
    p_z[5] = (sig_cur[5] * p_w[5] - sig_prv[5] * p_w[2]) / (-p_w[2] * p_w[2] - p_w[5] * p_w[5]);
    p_z[6] = (sig_cur[6] * p_w[6] - sig_prv[6] * p_w[1]) / (-p_w[1] * p_w[1] - p_w[6] * p_w[6]);
    p_z[7] = (sig_cur[7] * p_w[7] - sig_prv[7] * p_w[0]) / (-p_w[0] * p_w[0] - p_w[7] * p_w[7]);

    sig_prv[0] = p_work[3];
    sig_prv[1] = p_work[2];
    sig_prv[2] = p_work[1];
    sig_prv[3] = p_work[0];
    sig_prv[4] = -p_work[0];
    sig_prv[5] = -p_work[1];
    sig_prv[6] = -p_work[2];
    sig_prv[7] = -p_work[3];

サンプル数を可変

サンプル数を可変とすると下記になる。

    for (int i = 0; i < (nsmpl >> 1); i++) {
        sig_cur[i] = p_work[(nsmpl >> 1) + i];
        sig_cur[(nsmpl >> 1) + i] = p_work[nsmpl - 1 - i];
    }
    for (int i = 0; i < (nsmpl >> 1); i++) {
        p_z[i] = (sig_cur[i] * p_w[i] - sig_prv[i] * p_w[nsmpl - 1 - i]) /
                 (p_w[i] * p_w[i] + p_w[nsmpl - 1 - i] * p_w[nsmpl - 1 - i]);
        p_z[(nsmpl >> 1) + i] = -(sig_cur[(nsmpl >> 1) + i] * p_w[(nsmpl >> 1) + i] - sig_prv[(nsmpl >> 1) + i] * p_w[(nsmpl >> 1) - 1 - i]) /
                                 (p_w[(nsmpl >> 1) + i] * p_w[(nsmpl >> 1) + i] + p_w[(nsmpl >> 1) - 1 - i] * p_w[(nsmpl >> 1) - 1 - i]);
   }
   for (int i = 0; i < (nsmpl >> 1); i++) {
        sig_prv[i] = p_work[(nsmpl >> 1) - 1 - i];
        sig_prv[(nsmpl >> 1) + i] = -p_work[i];
   }

という感じです。

MDCTのコードを動かして動作を理解するのと、要素数を少なくして解析するということが重要でした。

一年前に良く分からなくて放置していた課題が解決してよかった。