I've run into an issue fiddling with the FFT/IFFT routines given here (http://www.dspguide.com/ch12/3.htm). I've rewritten the code to C++ and played with it. The idea is to read a raw wave file (encoded as 32 bit floats), find its FFT (complex) and then use that FFT as input to the InverseFFT function.
I expected to get the original sound data (samples) and I do for FFT lengths of up to 2048. For greater lengths (>= 4096 ), the code below doesn't (seem to) work. As you can see, my test consist solely of checking values at indexes 20,30 and 40 for both the original samples and the result (inverse FFT) array. I've tried this on several files (created with Audacity) with the same results.
int main(int argc, char *argv[])
{
unsigned fs;
float *Ix;
QCoreApplication a(argc, argv);
cout << "opening file...";
char *pData= OpenRaw("C:\\test\\250_350_c.raw",fs);
float* Rx = (float*)pData;
cout << "original data { [20,30,40] } " << Rx[20] << " " << Rx[30] << " " << Rx[40] << endl;
const unsigned numSam = 2048;
Ix=new float[numSam];
for(int i=0;i<numSam;i++)
Ix[i]= 0;
FFT(Rx,Ix,numSam);
float *Ifft = InverseFFT(Rx,Ix,numSam);
cout << "IFFT data { [20,30,40] } " << Ifft[20] << " " << Ifft[30] << " " << Ifft[40];
return a.exec();
}
char* OpenRaw(const char *fileName, unsigned& fs) {
char *pData;
unsigned fileSize = 0;
int bRead = 0;
ifstream file (fileName, ios::in | ios::binary);
if (file.is_open())
{
cout << "reading file..." << endl;
file.seekg (0, ios::end);
fileSize = file.tellg();
cout << "file size " << fileSize << endl;
pData = new char[fileSize];
file.seekg (0, ios::beg);
bRead = file.readsome(pData, fileSize);
if(file.fail())
cout << " error reading file..." << endl;
else
cout << "bytes read.." << bRead << endl;
file.close();
}else{
cout << "error reading file";
}
fs = fileSize / sizeof(float);
return pData;
}
float *InverseFFT(float *Rx, float *Ix, int fftSize){
float *pWaveData = new float[fftSize];
for(int i=0; i<fftSize; i++){
Ix[i] = -Ix[i];
}
FFT(Rx,Ix,fftSize);//calling the FFT function
for(int i=0; i<fftSize;i++){
pWaveData[i] = Rx[i]/fftSize;
Ix[i] = - Ix[i]/fftSize;
}
return pWaveData;
}
I haven't posted the FFT code here because I've used that code in the past and it always worked OK for me (I can still post it here if someone wishes to experiment with the code).
Does anyone know what might be causing this behaviour?
UPDATE 1: I've implemented the second method found in the link given by Peter below, but the result is still the same (the IFFT works for fft lengths of only up to 2048). My original IFFT method (given above) is actually the method #4 found in that same link.
float *InverseFFT(float *Rx, float *Ix, int fftSize){
FFT(Rx,Ix,fftSize);
float *pWaveData = new float[fftSize];
for(int i=1; i<fftSize/2; i++){
pWaveData[i] = Rx[fftSize-i]/fftSize;
pWaveData[fftSize-i] = Rx[i]/fftSize;
}
return pWaveData;
}
UPDATE 2: After googling some more, I found this ( https://www.mathworks.com/matlabcentral/newsreader/view_thread/247760?requestedDomain=www.mathworks.com ). Well, unfortunately, it seems that longer fft lengths may adversely influence an FFT/IFFT algorithm, especially if it's a power of 2 based algoritm.
Thank you