プログラミングメモ →目次

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

フィルターFIRの段数をNとして、一回で計算する入力信号の数をMとすると、N+M-1が2の整数乗になるようにNとMを決めます。 例えば、 FIRの段数 N=257、一回計算させる信号数 M = 768、FFTのサイズ = 1024、 重なる(Overlap)区域の信号数 = 1024 - 768 = 256 とします。

1024データ分の配列をゼロでクリアして、フィルターFIRを配列の左側にコピーし、FFTしておきます。入力信号をFFTする為に、1024データ分の配列を用意し、ゼロでクリアして、768個の信号データを配列にコピーします。入力信号を収納する配列をFFTし、周波数領域のフィルターと周波数領域の信号を複素乗算して、IFFTして結果を得ます。今回の結果(実数部分)の左端を、前回の結果の右端に、256データ分で重ねるようにして、加算(Add)して、最終の出力信号になります。

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

ソース:

#define FILTER_SAMPLE_NUM        257            //Sample numbers of LPF fir
#define SIGNAL_NUM_ALL           76800          //Sample numbers of signal to process
#define SAMPLE_FREQ              16000          //Sample frequency of signal to process
#define SIGNAL_NUM_PERSEGMENT    768            //Sample numbers per one segment
#define FFT_N                    1024           //Sample numbers for FFT
#define OVERLAP_NUM              256            //Sample numbers for overlap (= 1024 - 768)

extern double filter[FILTER_SAMPLE_NUM];     //LPF FIR Cutoff = 200Hz
double xr[FFT_N];                            //Array holding time domain signal real part for FFT                 
double xi[FFT_N];                            //Array holding time domain signal imaginary part for FFT                
double FR[FFT_N];                            //Array holding frequency domain filter real part for FFT                 
double FI[FFT_N];                            //Array holding frequency domain filter imaginary part for FFT                
double overlap[OVERLAP_NUM];                 //Array holding signal in overlap area 
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 OverlapAddMethodTest(CDC* pDC)
{
//Zero overlap array.
    ::ZeroMemory(overlap, sizeof(overlap));

//FFT for filter.
    ::ZeroMemory(FR, sizeof(FR));
    ::ZeroMemory(FI, sizeof(FI));
    ::CopyMemory(FR, filter, sizeof(filter));
    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; 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);
        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
        for (i = 0; i < FFT_N; i ++)
        {
            int x = SIGNAL_NUM_PERSEGMENT * (iOffset / SIGNAL_NUM_PERSEGMENT);
            if (i == 0)
                pDC->MoveTo((x+i)/2, 200+100 * (1 + iOffset / SIGNAL_NUM_PERSEGMENT)-(int)xr[i]);
            else
                pDC->LineTo((x+i)/2, 200+100 * (1 + iOffset / SIGNAL_NUM_PERSEGMENT)-(int)xr[i]);
        }

//Add overlap
        for (i = 0; i < OVERLAP_NUM; i ++)
        {
            xr[i] += overlap[i];
        }

//Prepare overlap for next

        for (i = SIGNAL_NUM_PERSEGMENT; i < FFT_N; i ++)
        {
            overlap[i - SIGNAL_NUM_PERSEGMENT] = xr[i];
        }
//Save output

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

    }
}



void COverlap_Add_MethodView::OnDraw(CDC* pDC)
{
    COverlap_Add_MethodDoc* 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 ++)
    {
        double alfa = 8.0 * atan(1.0) * 100 * i / SAMPLE_FREQ;
        input_signal[i] = 50.0 * sin(alfa) + 10.0 * sin(alfa * 8);  //Sine wave of 100hz with 800hz ripple.

    }
//Test.
    OverlapAddMethodTest(pDC);

//Display wave form.

    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, 200 - (int)output_signal[i]);
        else
            pDC->LineTo(i/2, 200  - (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
};