255個の係数のFIRフィルターのFFTによる実時間処理の実験

255個の係数のFIRフィルターのFFTによる実時間処理の実験

　SH-2A（144MHｚ　FPU付き）を使って、255個の係数をもつFIRフィルターの演算処理を実時間で実現できるかどうかを実験してみた。　オーディオ信号（サンプリング周波数44.1KHz 、ステレオ（２チャンネル）,16BIT整数、つまりCDと同じ）をSH-2Aの倍精度浮動小数点演算をつかい、演算量を減らして高速化するためオーバーラップアド法によるFFT演算を用い、SH-2Aの高速内臓RAMを用いて、PCからSH-2A基板へUSBハイスピードでオーディオデータを転送した場合、ほぼぎりぎりで実時間処理が可能なことが分かった。
但し、まだ、フィルターとしての減衰性能が十分に出ていない様子である。

■FFT OVERLAP ADD法
　
理論的に正しいかどうかは分からないが、演算量を少なくするため、
①ステレオ、２チャンネル信号を、一方のチャンネルを実数　もう片方のチャンネルを虚数として　FFTの複素数入力する。FFTからの出力も同様。　　FIRフィルターのインパルス応答が実数のみなので　つじつまがあうはずである？
②バタフライ演算をやらない。　正変換と逆変換の両方をおこなうので、入れ替え（バタフライ演算）はしなくてももとの順番に戻せるはず？
としている。

叉、上記①の計算をしたためか、単精度浮動小数点演算では、処理時間はおよそ半分ぐらいになるが、音が（両チャンネルの干渉が大きくでるためか？）良くないので、採用しなかった。

FFTのCOS,SINテーブルや頻繁に使う作業用の記憶域が、SH-2Aの高速RAM上にないと、アクセルに時間がかかるため、FFTの計算時間が長くなり、実時間ではできない。しかし、SH-2Aの高速RAMの容量サイズ制限のため、すべてのデータはそこにはおけない（すでのSTACK領域も使いこんでいる）。そこで、残ったFIRフィルターをFFTしたデータはSH-2Aのキャッシュ（８K)においてもらうように、唯一このデータだけ、大容量 RAMのキャッシュ対応アドレスに割り当てた（上手く行っているかどうかは分からないが・・・）

// 全データ領域をを高速RAMのセッションにすると起動しなくなる。
// STACK領域を使い込んでいるので注意。
#pragma section WDATA //このセッションを高速RAMに定義する
double mk[klen0 * 2];    // cos,sinの係数テーブル        //8k
double works[klen0 * 2]; // 作業領域       //8k
//　ここまででスタック領域に進入している　stack の残り約12k
#pragma section

double works_half[klen0_half * 2];

#pragma section CDATA //大容量RAMのキャッシュ対応アドレスのセッションとして定義する。
double kaiser_lpf2[klen0 * 2];   // FIR フィルターの応答をFFTしたものを格納しておく。
#pragma section

int f_fft_overlap_add2(short *in0, int count0)
{
   int n,nh,i,i1,i2,x1,x2,offset;
   double *ws,*datas,a0,a1;
   short *in0p;

   /*　16BIT 整数信号をDOBLUEへ変換 */
if( count0 >= 0)
{
   datas=works;
   memcpy((char *) datas, (char *) works_half, klen0 * sizeof(double) );
   in0p=in0;
   datas=works+klen0;
   for(i=klen0_half;i<klen0;++i)
   {
        *datas   = (double) *in0p;
        ++datas;
        ++in0p;
        *datas   = (double) *in0p;
        ++datas;
        ++in0p;
   }
   // 次回のためにデータを退避

    datas=works;
    memcpy((char *) works_half, ((char *)datas + klen0 * sizeof(double)) , klen0 * sizeof(double));

} //if( count0 >= 0)
else
// FIR の応答のFFTを計算する。
{
    datas=kaiser_lpf2;
}

   /* 正変換 */
   /* start fft */
   nh=klen0_half;
   n=1;
   ws=mk;
   for(i1=0;i1<nzyoux;++i1)
   {
         for(i2=0;i2<n;++i2)     // f_fttsub()で2個分　計算しているので　1個少ない。
         {

            offset=2 * nh*i2;   // nh　ステップのi2番目
            for(i=0;i<nh;++i)
            {
               x1=2*(i+offset);
               x2=2*(i+nh+offset);
               a0= datas[x1] - datas[x2];
               a1= datas[x1+1] - datas[x2+1];
               datas[x1]   = datas[x1] + datas[x2];
               datas[x1+1] = datas[x1+1] + datas[x2+1];

               // 2個計算している
               datas[x2]   = a0 * ws[2*i*n]   - a1 * ws[2*i*n+1];
               datas[x2+1] = a0 * ws[2*i*n+1] + a1 * ws[2*i*n];
             }
         }
         nh=nh>>1; //nh=nh/2;
         n=n<<1;    //n=n*2;
   }

/* 正変換のみの場合 */
if( count0 < 0) return(999);

   /* 周波数領域で乗算する */
for(i=0;i<klen0;++i)
{
     a0= datas[2*i];
      datas[2*i]   = datas[2*i] * kaiser_lpf2[2*i] - datas[2*i+1] * kaiser_lpf2[2*i+1];
     datas[2*i+1] = datas[2*i+1] * kaiser_lpf2[2*i] + a0 * kaiser_lpf2[2*i+1];
}

/* 逆変換 */
/* start fft */
// wk の虚数部を符号逆にしてつかう。

nh=1;
n=klen0_half;
for(i1=0;i1<nzyoux;++i1)
{
         for(i2=0;i2<n;++i2)
         {

            offset=2*nh*i2;
            for(i=0;i<nh;++i)
            {
               x1=2*(i+offset);
               x2=2*(i+nh+offset);
               a0   = datas[x2] * ws[2*i*n]   + datas[x2+1] * ws[2*i*n+1];
               a1 = - datas[x2] * ws[2*i*n+1] + datas[x2+1] * ws[2*i*n];
               datas[x2]   = datas[x1] - a0;
               datas[x2+1] = datas[x1+1] - a1;
               datas[x1]   = datas[x1] + a0;
               datas[x1+1] = datas[x1+1] + a1;
            }

         }
         n=n>>1;    //n=n/2;
         nh=nh<<1; //nh=nh * 2;
    }

   /* 1/N */

   in0p=in0;
   datas= works+klen0;
   for(i=(klen0_half);i<klen0;++i)   // 前半は捨てる
   {
        // a0 = datas[2*i];
        a0= *datas;
        if(a0 >=0)
        {
             if(a0 > 32767.) a0=32767.;
             *in0p= (short)(a0 + 0.5);   // 4捨5入
        }
        else
        {
             if(a0 < -32768.) a0=-32768.;
             *in0p= (short)(a0 - 0.5); // 4捨5入
        }
        ++datas;
        ++in0p;
        // a0 = datas[2*i+1];
        a0= *datas;
        if(a0 >=0)
        {
             if(a0 > 32767.) a0=32767.;
             *in0p= (short)(a0 + 0.5);   // 4捨5入
        }
        else
        {
             if(a0 < -32768.) a0=-32768.;
             *in0p= (short)(a0 - 0.5); // 4捨5入
        }
        ++datas;
        ++in0p;
   }
   return(0);
}

■ハード
　先回、PCからSH-2A基板へオーディオデータをバルク転送する実験をおこなった基板に、更に、SH-2Aの中でFFT OVERLAP ADD法の計算をしている時間が分かるように、SH-2AのPG14ポート(CN3の17番)に黄色（橙色）のLEDを光らせる回路（先回つけた赤LEDと同様回路）をつけて、FFT　OVERLAP　 ADD法を計算している間、LEDを点灯させることにした。

参考用３: main.zip 上記のソースファイル　 PC側のソフトから一度STOP　BITを０に設定した後にSTOP　BITを２に設定すると　FFT　OVERLAP　ADD法の計算をするようになっている。FIRフィルターには、255個の係数をもつ、仮のフィルター（ほとんど通過域のローパスフィルター）の値が入っている。

（追加）FFT OVERLAP ADD法による4倍のオーバーサンプリング用デジタルフィルタを作ってみた。

サンプリングされた入力の間に零を３個入れた信号に、デジタルフィルターを掛けることにより原信号の波形の形に近づけようとするものである（完全には同じにならない）。
512個の係数をもつFIRフィルターを、零を内挿したため、128(129)個の係数をもつ4組のFIRフィルターに分割できる。カットオフ周波数をサンプリング周波数（この場合は、44100Hzの4倍の176400Hz）の1/8に選ぶと、はじめの組は中心以外の係数が零になるため、入力信号を64(256)サンプル　シフトさせたものになり、計算が省略できる。零を内挿して4個あるべきところ1個しかない歯抜け波形になって、元の信号波形に対してパワーが４分の１になってしまったので、残りの3組のFIRフィルターの出力には４倍を掛けている。
カットオフ周波数での減衰量は、理論値で約－6ｄBである。減衰目標値は－110ｄBであるが、実際には当てにならない。本当は、倍精度で計算したかったのであるが、実時間処理できなかったので、仕方なく、単精度の浮動小数点演算で計算している。180度位相がずれた１KHｚのSIN波入力して計算させた結果、（16ビットですでに）下位の1桁に誤差が見られた。なので、演算精度は15ビット程度以下しかないようだ。全体的に今一の出来のフィルターとなっている。なを、参考用に内容を上記のmain.zipに入れてある。PC側のソフトから一度STOP　BITを０に設定した後にSTOP　BITを４に設定する。4倍のオーバーサンプリング用デジタルフィルタするようになっている。

■ 警告
SH-2A基板をつかって接続する場合は、電気的なことをよく理解した上でお使いください。そうしないと、PCやその基板と接続した相手などを破壊する危険があります。
供給元が違うVBUS＋５V同士を接続すると、電源が喧嘩をして破壊する恐れがあります。供給元が違う電源は一緒に接続しないように。

■免責
(1)回路図やプログラムやデータの使用により、使用者に損失が生じたとしても、その責任を負いません。
(2)回路図やプログラムやデータにバグや欠陥があったとしても、修正や改良の義務を負いません。