求频率细化分析zoomfft的C++源代码

求频率细化分析zoomfft的C++源代码

//下面的FFT我用了很多年了：

// 离散傅里叶变换DFT代码：
int DFT (long count, CComplex * input, CComplex * output)
{
assert(count);
assert(input);
assert(output);

CComplex F, X, T, W; int n, i;

long N = abs(count); long Inversing = count < 0? 1: -1;

for(n = 0; n < N ; n++){ // compute from line 0 to N-1

F = CComplex(0.0f, 0.0f); // clear a line

for(i = 0; i < N; i++) {

T = input[i];

W = HarmonicPI2(Inversing * n * i, N);

X = T * W;

F += X; // fininshing a line

}//next i

// save data to outpus
memcpy(output + n, &F, sizeof(F));

}//next n

return 0;
}//end DFT

//快速傅里叶变换代码FFT
int fft (long count, CComplex * input, CComplex * output)
{
assert(count);
assert(input);
assert(output);

int N = abs(count); long Inversing = count < 0? -1: 1;

if (N % 2 || N < 5) return DFT(count, input, output);

long N2 = N / 2;

CComplex * iEven = new CComplex[N2]; memset(iEven, 0, sizeof(CComplex) * N2);
CComplex * oEven = new CComplex[N2]; memset(oEven, 0, sizeof(CComplex) * N2);
CComplex * iOdd = new CComplex[N2]; memset(iOdd , 0, sizeof(CComplex) * N2);
CComplex * oOdd = new CComplex[N2]; memset(oOdd , 0, sizeof(CComplex) * N2);

int i = 0; CComplex W;
for(i = 0; i < N2; i++) {
iEven[i] = input[i * 2];
iOdd [i] = input[i * 2 + 1];
}//next i

fft(N2 * Inversing, iEven, oEven);
fft(N2 * Inversing, iOdd, oOdd );

for(i = 0; i < N2; i++) {
W = HarmonicPI2(Inversing * (- i), N);
output[i] = oEven[i] + W * oOdd[i];
output[i + N2] = oEven[i] - W * oOdd[i];
}//next i
return 0;
}//end FFT