FFT畳み込み積分(FFT Convolution), Overlap-Add Method フィルターFIRの段数をNとして、一回で計算する入力信号の数をMとすると、N+M-1が２の整数乗になるように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
};
|
