FFT畳み込み積分(FFT Convolution), Overlap-Save Method

FIRフィルターの段数をNとして、一回に計算する入力信号数を２Mとすると、M >= N 　且つ 　２Mが２の整数乗になるようにMを決めます。 例えば、 FIRの段数 N=257、一回計算させる入力信号数 M = 512、FFTのサイズ = 1024　とします。フィルターを収納する配列のサイズを２Mにして、中心から左側の値をゼロにして、中心から右側にFIRフィルターの係数を代入します。右端に余った部分にもゼロを代入し、そしてFFTします。１０２４個の入力信号をFFTします。周波数領域のフィルターと周波数領域の入力信号を複素乗算して、逆FFTして時間領域信号を求めます。逆FFT の結果の時間領域信号の前半半分が直線状畳み込み結果となるので、前半半分をつなげていくことで畳み込みを実現します。時間領域信号の後半半分は、「円状状畳み込み結果」なので、使いません。

入力信号をM個進めて、次の計算を行います。同じように、M個の出力信号を得ます。隣り合いのIFFTで得た時間領域信号の前半半分をつなげていくことで畳み込みを実現します。

例：サンプリング周波数の16000HZを持つ入力信号に、Cutoff周波数の200HZのLPFフィルターを掛けます。入力信号には、100HZの正弦波と小さい800HZ正弦波が混じっています。

ソース：

#define FILTER_SAMPLE_NUM       257     //Sample numbers of LPF fir
#define SIGNAL_NUM_ALL          25600   //Sample numbers of signal to process
#define SAMPLE_FREQ             16000   //Sample frequency of signal to process
#define SIGNAL_NUM_PERSEGMENT   512     //Sample numbers per one segment
#define FFT_N                   1024    //Sample numbers for FFT

extern double filter[FILTER_SAMPLE_NUM]; //LPF FIR Cutoff = 200Hz
double xr[FFT_N];                        //Array holding time domain signal real part            
double xi[FFT_N];                        //Array holding time domain signal imaginary part                
double FR[FFT_N];                        //Array holding frequency domain filter real part                 
double FI[FFT_N];                        //Array holding frequency domain filter imaginary part                
double* input_signal;                    //Array holding input signal                
double* output_signal;                   //Array holding output signal                

void fft(int n, double theta, double ar[], double ai[]);
void OverlapSaveMethodTest(CDC* pDC)
{
    CPen penGray(PS_SOLID, 1, (COLORREF)0x00c0c0c0);
    CPen penBlack(PS_SOLID, 1, (COLORREF)0x00000000);

//FFT for filter.
    ::ZeroMemory(FR, sizeof(FR));
    ::ZeroMemory(FI, sizeof(FI));
    ::CopyMemory(&FR[FFT_N/2], filter, sizeof(double)*min(FFT_N/2, FILTER_SAMPLE_NUM));
    double theta = -8 * atan(1.0) / FFT_N;
    fft(FFT_N, theta, FR, FI);
//Loop for each segment
    for (int iOffset = 0; iOffset < SIGNAL_NUM_ALL - SIGNAL_NUM_PERSEGMENT; 
                iOffset += SIGNAL_NUM_PERSEGMENT)
    {
//FFT for one segment
        ::ZeroMemory(xr, sizeof(xr));
        ::ZeroMemory(xi, sizeof(xi));
        ::CopyMemory(xr, &input_signal[iOffset], sizeof(double) * SIGNAL_NUM_PERSEGMENT * 2);
        double theta = -8 * atan(1.0) / FFT_N;
        fft(FFT_N, theta, xr, xi);
//Multiply frequency spectrum of filter and frequency spectrum of signal
        for (int i = 0; i < FFT_N; i ++)
        {
            double dbl_R = FR[i] * xr[i] - FI[i] * xi[i];
            double dbl_I = FR[i] * xi[i] + FI[i] * xr[i];
            xr[i] = dbl_R;
            xi[i] = dbl_I;
        }
//IFFT for one segment
        theta = 8 * atan(1.0) / FFT_N;
        fft(FFT_N, theta, xr, xi);
        for (i = 0; i < FFT_N; i ++)
        {
            xr[i] /= FFT_N;
        }

//Display xr
        pDC->SelectObject(&penBlack);
        for (i = 0; i < FFT_N; i ++)
        {
            if (i >= SIGNAL_NUM_PERSEGMENT)
            {
                pDC->SelectObject(&penGray);
            }
            int x = SIGNAL_NUM_PERSEGMENT * (iOffset / SIGNAL_NUM_PERSEGMENT);
            if (i == 0)
                pDC->MoveTo((x+i)/2, 100+100 * (1 + iOffset / SIGNAL_NUM_PERSEGMENT)-(int)xr[i]);
            else
                pDC->LineTo((x+i)/2, 100+100 * (1 + iOffset / SIGNAL_NUM_PERSEGMENT)-(int)xr[i]);
        }

//Save output

        ::CopyMemory(&output_signal[iOffset], xr, sizeof(double) * SIGNAL_NUM_PERSEGMENT);

    }
}

void COverlapSaveMethodView::OnDraw(CDC* pDC)
{
    COverlapSaveMethodDoc* pDoc = GetDocument();
    ASSERT_VALID(pDoc);
    // TODO: この場所にネイティブ データ用の描画コードを追加します。

//Make signals to test.
    input_signal = new double [SIGNAL_NUM_ALL];
    output_signal = new double [SIGNAL_NUM_ALL];
    double pi2 = 8.0 * atan(1.0);
    for (int i = 0; i < SIGNAL_NUM_ALL; i ++)
    {
        //Sine wave of 100hz with 800hz ripple.
        double alfa = 8.0 * atan(1.0) * 100 * i / SAMPLE_FREQ;
        input_signal[i] = 50.0 * sin(alfa) + 10.0 * sin(alfa * 8);
    }
//Test.
    OverlapSaveMethodTest(pDC);

//Display wave form.
    CPen penBlack(PS_SOLID, 1, (COLORREF)0x00000000);
    pDC->SelectObject(&penBlack);
    for (i = 0; i < 3200; i ++)
    {
        if (i == 0)
            pDC->MoveTo(i/2, 100 - (int)input_signal[i]);
        else
            pDC->LineTo(i/2, 100  - (int)input_signal[i]);
    }

    for (i = 0; i < 3200; i ++)
    {
        if (i == 0)
            pDC->MoveTo(i/2, 600 - (int)output_signal[i]);
        else
            pDC->LineTo(i/2, 600  - (int)output_signal[i]);
    }

    delete [] input_signal;
    delete [] output_signal;
}






void fft(int n, double theta, double ar[], double ai[])
{
    int m, mh, i, j, k;
    double wr, wi, xr, xi;

    for (m = n; (mh = m >> 1) >= 1; m = mh) 
    {
        for (i = 0; i < mh; i++) 
        {
            wr = cos(theta * i);
            wi = sin(theta * i);
            for (j = i; j < n; j += m) 
            {
                k = j + mh;
                xr = ar[j] - ar[k];
                xi = ai[j] - ai[k];
                ar[j] += ar[k];
                ai[j] += ai[k];
                ar[k] = wr * xr - wi * xi;
                ai[k] = wr * xi + wi * xr;
            }
        }
        theta *= 2;
    }
    /* ---- unscrambler ---- */
    i = 0;
    for (j = 1; j < n - 1; j++) 
    {
        for (k = n >> 1; k > (i ^= k); k >>= 1);
        if (j < i) 
        {
            xr = ar[j];
            xi = ai[j];
            ar[j] = ar[i];
            ai[j] = ai[i];
            ar[i] = xr;
            ai[i] = xi;
        }
    }
}

double filter[FILTER_SAMPLE_NUM] =
{
-0.00019086,
-0.000195,
-0.00019961,
-0.00020469,
-0.00021027,
-0.00021636,
-0.00022296,
-0.00023008,
-0.00023771,
-0.00024583,
-0.00025441,
-0.00026344,
-0.00027287,
-0.00028265,
-0.00029274,
-0.00030308,
-0.00031359,
-0.00032421,
-0.00033485,
-0.00034541,
-0.0003558,
-0.00036592,
-0.00037565,
-0.00038487,
-0.00039346,
-0.00040127,
-0.00040818,
-0.00041404,
-0.0004187,
-0.000422,
-0.00042379,
-0.00042389,
-0.00042215,
-0.00041839,
-0.00041243,
-0.0004041,
-0.00039323,
-0.00037963,
-0.00036313,
-0.00034355,
-0.00032071,
-0.00029443,
-0.00026454,
-0.00023086,
-0.00019324,
-0.0001515,
-0.00010548,
-5.5028e-005,
1.789e-019,
5.975e-005,
0.00012436,
0.00019395,
0.00026864,
0.00034855,
0.00043376,
0.00052439,
0.00062049,
0.00072215,
0.00082941,
0.00094233,
0.0010609,
0.0011852,
0.0013153,
0.001451,
0.0015924,
0.0017395,
0.0018921,
0.0020503,
0.002214,
0.002383,
0.0025573,
0.0027367,
0.0029212,
0.0031104,
0.0033044,
0.0035028,
0.0037056,
0.0039125,
0.0041232,
0.0043376,
0.0045554,
0.0047763,
0.005,
0.0052264,
0.005455,
0.0056856,
0.0059179,
0.0061515,
0.0063861,
0.0066214,
0.006857,
0.0070926,
0.0073278,
0.0075622,
0.0077955,
0.0080273,
0.0082572,
0.0084849,
0.0087099,
0.008932,
0.0091507,
0.0093657,
0.0095765,
0.0097829,
0.0099845,
0.010181,
0.010372,
0.010557,
0.010736,
0.010908,
0.011073,
0.011232,
0.011383,
0.011526,
0.011661,
0.011788,
0.011907,
0.012017,
0.012118,
0.01221,
0.012293,
0.012367,
0.012431,
0.012485,
0.012529,
0.012564,
0.012589,
0.012604,
0.012609,
0.012604,
0.012589,
0.012564,
0.012529,
0.012485,
0.012431,
0.012367,
0.012293,
0.01221,
0.012118,
0.012017,
0.011907,
0.011788,
0.011661,
0.011526,
0.011383,
0.011232,
0.011073,
0.010908,
0.010736,
0.010557,
0.010372,
0.010181,
0.0099845,
0.0097829,
0.0095765,
0.0093657,
0.0091507,
0.008932,
0.0087099,
0.0084849,
0.0082572,
0.0080273,
0.0077955,
0.0075622,
0.0073278,
0.0070926,
0.006857,
0.0066214,
0.0063861,
0.0061515,
0.0059179,
0.0056856,
0.005455,
0.0052264,
0.005,
0.0047763,
0.0045554,
0.0043376,
0.0041232,
0.0039125,
0.0037056,
0.0035028,
0.0033044,
0.0031104,
0.0029212,
0.0027367,
0.0025573,
0.002383,
0.002214,
0.0020503,
0.0018921,
0.0017395,
0.0015924,
0.001451,
0.0013153,
0.0011852,
0.0010609,
0.00094233,
0.00082941,
0.00072215,
0.00062049,
0.00052439,
0.00043376,
0.00034855,
0.00026864,
0.00019395,
0.00012436,
5.975e-005,
1.789e-019,
-5.5028e-005,
-0.00010548,
-0.0001515,
-0.00019324,
-0.00023086,
-0.00026454,
-0.00029443,
-0.00032071,
-0.00034355,
-0.00036313,
-0.00037963,
-0.00039323,
-0.0004041,
-0.00041243,
-0.00041839,
-0.00042215,
-0.00042389,
-0.00042379,
-0.000422,
-0.0004187,
-0.00041404,
-0.00040818,
-0.00040127,
-0.00039346,
-0.00038487,
-0.00037565,
-0.00036592,
-0.0003558,
-0.00034541,
-0.00033485,
-0.00032421,
-0.00031359,
-0.00030308,
-0.00029274,
-0.00028265,
-0.00027287,
-0.00026344,
-0.00025441,
-0.00024583,
-0.00023771,
-0.00023008,
-0.00022296,
-0.00021636,
-0.00021027,
-0.00020469,
-0.00019961,
-0.000195,
-0.00019086
};