📄 fft_radix4_2_cplx.c
字号:
/*
* DATE:05/29/2006
*
* AUTHOR:CHEN PENG
*
* VERSION:ORIGINAL 05/29/2006
* MODIFIED 06/07/2006
* 06/23/2006
*/
/*****************************************************************************************\
* FUNCTION: void fft_radix4_2_cplx( int x[], int num_r2, int w[], short power0f4)
*
* AGRUMENTS: x[num_r2] Input and output sequences. Size num_r2 complex
* elements
*
* num_r2 Number of complex elements in vector x. Must be
* power of 2 such that 4 <= num_r2 <= 65536
*
* w[num_r2] Vector of FFT coefficients of size num_r2 elements
*
* power0f4 Flag which indicates if size of sequence is power of 4
*
* DESCRIPTION: This routine is used to compute FFT of a complex sequence
* size num_r2, a power of 2 or 4, with "decimation-in-frequency
* radix4 combined with radix2 decomposition" method . The
* output is in digit-reversed order. Each complex value is with
* interleaved 16-bit real and imaginary parts.
*
* REQUIREMENTS: 1. 4 <= num_r2 <= 65536( num_r2 a power of 2)
* 2. x data is stored in the order {image[0],real[0]},{image[1],...
* 3. x is aligned on a 4*num_r2 byte boundary
* 4. The code is LITTLE ENDIAN
*
\*****************************************************************************************/
void fft_radix4_2_cplx( int x[], int num_r2, int w[], short powerOf4 ){
int n1, n2, ie, num_r4; //Loop counters
int iw1, iw2, iw3; //FFT coefficients index
int ix0, ix1, ix2, ix3; //Data sequences index
int ix0b, ix1b, ix2b, ix3b; //Data sequences index only for radix4 mixed with radix2
int iStage, iButterfly; //Stage and butterfly flag
int xtmph, xtmpl; //Temporary high and low 16-bit results
int temp00, temp01, temp02; //Temporary variables
int temp10, temp11, temp12;
num_r4 = num_r2;
n2 = num_r2;
ie = 1;
/*
*If input sequences has size of power of 4, FFT will be computed with method radix4,
*othewise, firstly the input sequences is decomposited with method radix2 into two
*separate sequences each with size of num_r4, then FFT calculation for each part
*performs under radix4 respectively.
*/
if( !powerOf4 ){
/*===================*\
*radix2 FFT firstly
\*===================*/
num_r4 = num_r2 >> 1; //size of each part for radix4 method
n1 = n2;
n2 >>= 1;
iw1 = 0;
for(iButterfly = 0; iButterfly < n2; iButterfly++){//number of butterflies
for(ix0 = iButterfly; ix0 < num_r2; ix0 += n1){//loop for butterfly calculations
ix2 = ix0 + n2;
temp00 = _add2(x[ix0],x[ix2]);
temp01 = _sub2(x[ix0],x[ix2]);
x[ix0] = temp00; //even points
xtmph = (_smpyh(temp01,w[iw1]) - _smpy(temp01,w[iw1])) & 0xffff0000;
xtmpl = ((_smpylh(temp01,w[iw1]) + _smpyhl(temp01,w[iw1])) >> 16)
& 0x0000ffff;
x[ix2] = xtmph | xtmpl; //odd points
}
iw1 += 1;
}
/*======================*\
*radix4 FFT afterwards
\*======================*/
/*
* X(4r) = FFT[A(n)] A(n) = x(n) + x(n + N/2) + x(n + N/4) + x(n + 3N/4)
* X(4r + 2) = FFT[B(n)] B(n) = [x(n) + x(n + N/2) - x(n + N/4) - x(n + 3N/4)]*Wn(2n)
* X(4r + 1) = FFT[C(n)] C(n) = [x(n) - x(n + N/2) - jx(n + N/4) + jx(n + 3N/4)]*Wn(n)
* X(4r + 3) = FFT[D(n)] D(n) = [x(n) - x(n + N/2) + jx(n + N/4) - jx(n + 3N/4)]*Wn(3n)
*/
for( iStage = num_r4; iStage > 1; iStage >>= 2){//number of stages
n1 = n2;
n2 >>= 2;
iw1 = 0;
for(iButterfly = 0; iButterfly < n2; iButterfly++){//number of butterflies
iw2 = iw1 + iw1;
iw3 = iw2 + iw1;
for(ix0 = iButterfly; ix0 < num_r4; ix0 += n1){//loop for butterfly calculations
ix1 = ix0 + n2;
ix2 = ix1 + n2;
ix3 = ix2 + n2;
/****************************************************************\
*even points calculations for former sequence with num_r4 points
\****************************************************************/
temp00 = _add2(x[ix1],x[ix3]);
temp01 = _add2(x[ix0],x[ix2]);
temp02 = _sub2(x[ix0],x[ix2]);
x[ix0] = _add2(temp00,temp01); //A(n)
temp01 = _sub2(temp01,temp00);
xtmph = (_smpyh(temp01,w[iw2]) - _smpy(temp01,w[iw2])) & 0xffff0000;
xtmpl = ((_smpylh(temp01,w[iw2]) + _smpyhl(temp01,w[iw2])) >> 16)
& 0x0000ffff;
temp00 = _sub2(x[ix1],x[ix3]);
x[ix1] = xtmph | xtmpl; //B(n)
/****************************************************************\
*odd points calculations for former sequence with num_r4 points
\****************************************************************/
temp01 = (temp00 << 16);
temp00 = temp01 | ((-(temp00 >> 16)) & 0x0000ffff);
temp01 = _add2(temp02,temp00);
temp02 = _sub2(temp02,temp00);
xtmph = (_smpyh(temp01,w[iw1]) - _smpy(temp01,w[iw1])) & 0xffff0000;
xtmpl = ((_smpylh(temp01,w[iw1]) + _smpyhl(temp01,w[iw1])) >> 16)
& 0x0000ffff;
x[ix2] = xtmph | xtmpl; //C(n)
xtmph = (_smpyh(temp02,w[iw3]) - _smpy(temp02,w[iw3])) & 0xffff0000;
xtmpl = ((_smpylh(temp02,w[iw3]) + _smpyhl(temp02,w[iw3])) >> 16)
& 0x0000ffff;
x[ix3] = xtmph | xtmpl; //D(n)
/****************************************************************\
*even points calculations for latter sequence with num_r4 points
\****************************************************************/
ix0b = ix0 + num_r4;
ix1b = ix1 + num_r4;
ix2b = ix2 + num_r4;
ix3b = ix3 + num_r4;
temp10 = _add2(x[ix1b],x[ix3b]);
temp11 = _add2(x[ix0b],x[ix2b]);
temp12 = _sub2(x[ix0b],x[ix2b]);
x[ix0b] = _add2(temp10,temp11); //A(n)
temp11 = _sub2(temp11,temp10);
xtmph = (_smpyh(temp11,w[iw2]) - _smpy(temp11,w[iw2])) & 0xffff0000;
xtmpl = ((_smpylh(temp11,w[iw2]) + _smpyhl(temp11,w[iw2])) >> 16)
& 0x0000ffff;
temp10 = _sub2(x[ix1b],x[ix3b]);
x[ix1b] = xtmph | xtmpl; //B(n)
/****************************************************************\
*odd points calculations for latter sequence with num_r4 points
\****************************************************************/
temp11 = (temp10 << 16);
temp10 = temp11 | ((-(temp10 >> 16)) & 0x0000ffff);
temp11 = _add2(temp12,temp10);
temp12 = _sub2(temp12,temp10);
xtmph = (_smpyh(temp11,w[iw1]) - _smpy(temp11,w[iw1])) & 0xffff0000;
xtmpl = ((_smpylh(temp11,w[iw1]) + _smpyhl(temp11,w[iw1])) >> 16)
& 0x0000ffff;
x[ix2b] = xtmph | xtmpl; //C(n)
xtmph = (_smpyh(temp12,w[iw3]) - _smpy(temp12,w[iw3])) & 0xffff0000;
xtmpl = ((_smpylh(temp12,w[iw3]) + _smpyhl(temp12,w[iw3])) >> 16)
& 0x0000ffff;
x[ix3b] = xtmph | xtmpl; //D(n)
}
iw1 += (2 * ie); //Scaler '2' is used when radix2 computation
} //is performed first
ie <<= 2;
}
}
else{
/*==================*\
*direct radix4 FFT
\*==================*/
for( iStage = num_r4; iStage > 1; iStage >>= 2){
n1 = n2;
n2 >>= 2;
iw1 = 0;
for(iButterfly = 0; iButterfly < n2; iButterfly++){
iw2 = iw1 + iw1;
iw3 = iw2 + iw1;
for(ix0 = iButterfly; ix0 < num_r4; ix0 += n1){
ix1 = ix0 + n2;
ix2 = ix1 + n2;
ix3 = ix2 + n2;
/*************************\
*even points calculations
\*************************/
temp00 = _add2(x[ix1],x[ix3]);
temp01 = _add2(x[ix0],x[ix2]);
temp02 = _sub2(x[ix0],x[ix2]);
x[ix0] = _add2(temp00,temp01); //A(n)
temp01 = _sub2(temp01,temp00);
xtmph = (_smpyh(temp01,w[iw2]) - _smpy(temp01,w[iw2])) & 0xffff0000;
xtmpl = ((_smpylh(temp01,w[iw2]) + _smpyhl(temp01,w[iw2])) >> 16)
& 0x0000ffff;
temp00 = _sub2(x[ix1],x[ix3]);
x[ix1] = xtmph | xtmpl; //B(n)
/*************************\
*odd points calculations
\*************************/
temp01 = (temp00 << 16);
temp00 = temp01 | ((-(temp00 >> 16)) & 0x0000ffff);
temp01 = _add2(temp02,temp00);
temp02 = _sub2(temp02,temp00);
xtmph = (_smpyh(temp01,w[iw1]) - _smpy(temp01,w[iw1])) & 0xffff0000;
xtmpl = ((_smpylh(temp01,w[iw1]) + _smpyhl(temp01,w[iw1])) >> 16)
& 0x0000ffff;
x[ix2] = xtmph | xtmpl; //C(n)
xtmph = (_smpyh(temp02,w[iw3]) - _smpy(temp02,w[iw3])) & 0xffff0000;
xtmpl = ((_smpylh(temp02,w[iw3]) + _smpyhl(temp02,w[iw3])) >> 16)
& 0x0000ffff;
x[ix3] = xtmph | xtmpl; //D(n)
}
iw1 += ie;
}
ie <<= 2;
}
}
return;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -