Problem getting IFFT to work

Question

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

Answer 1

I solved this by using a different FFT implementation, one that is using a mixed radix algorithm under the hood ( the project is avalaible here https://github.com/itdaniher/kissfft ). The KissFFT supports both a complex as well as a real FFT so I'll post code for both. The real version is up to 45% faster than the complex one. Note that in order to get the original samples, the results of the IFFT have to be scaled by fft size (fft size doesn't have to be a power of 2 but for the real FFT it must be an even number).

The complex FFT/IFFT

    #include "kiss_fft.h"
    #include "kiss_fftr.h"

        int main(int argc, char *argv[])
            {

                unsigned fs;

                QCoreApplication a(argc, argv);
                cout << "opening ...";
                char *pData= OpenRawA("C:\\test\\250_350_c.raw",fs);

                const unsigned numSamples = 4096;//this doesn't have to be a power of 2
                float* Rx = (float*)pData;


                       kiss_fft_cpx *fIn = new kiss_fft_cpx[numSamples];
                       kiss_fft_cpx *fOut = new kiss_fft_cpx[numSamples] ;
                       kiss_fft_cpx *fIfft = new kiss_fft_cpx[numSamples];
                       //populate the (complex) input buffer
                        for(int i=0; i<numSamples;i++){
                            fIn[i].r = Rx[i];
                            fIn[i].i = 0;
                        }


                        cout << "original data { [20,30,40] } " << fIn[20].r << " " << fIn[30].r << " " << fIn[40].r <<  fIn[50].r << " " << fIn[60].r <<endl;
                        kiss_fft_cfg cfg = kiss_fft_alloc( numSamples ,0 ,0,0 );//configure the FFT function to use a forward FFT

                        kiss_fft( cfg , fIn , fOut );//calculate FFT


                        cfg = kiss_fft_alloc( numSamples ,1 ,0,0 );//configure the FFT function to use an inverse FFT

                        kiss_fft( cfg , fOut , fIfft );//calculate inverse FFT

                       free(cfg);


                cout << "IFFT data { [20,30,40] } " << fIfft[20].r/numSamples << " " << fIfft[30].r/numSamples << " " << fIfft[40].r/numSamples << " " << fIfft[50].r/numSamples << " " << fIfft[60].r/numSamples;
//delete allocated resources
                return a.exec();
            }

The real FFT/IFFT

 int main(int argc, char *argv[])
    {

        unsigned fs;

        QCoreApplication a(argc, argv);
        cout << "opening ...";
         char *pData= OpenRawA("C:\\test\\250_350_c.raw",fs);

        const unsigned numSamples = 4096;
        float* Rx = (float*)pData;

               kiss_fft_scalar *fIn = new kiss_fft_scalar[numSamples];
               kiss_fft_cpx *fOut = new kiss_fft_cpx[numSamples] ;
               kiss_fft_scalar *fIfft = new kiss_fft_scalar[numSamples];
               //populate the input buffer
                for(int i=0; i<numSamples;i++){
                    fIn[i] = Rx[i];

                }
                cout << "original data { [20,30,40] } " << fIn[20] << " " << fIn[30] << " " << fIn[40] <<  fIn[50] << " " << fIn[60] <<endl;
                kiss_fftr_cfg cfg = kiss_fftr_alloc( numSamples ,0 ,0,0 );//configure the FFT function to use a forward FFT

                kiss_fftr( cfg , fIn , fOut );//calculate FFT

                cfg = kiss_fftr_alloc( numSamples ,1 ,0,0 );//configure the FFT function to use an inverse FFT

                kiss_fftri( cfg , fOut , fIfft );//calculate inverse FFT

               free(cfg);

        cout << "IFFT data { [20,30,40] } " << fIfft[20]/numSamples << " " << fIfft[30]/numSamples << " " << fIfft[40]/numSamples << " " << fIfft[50]/numSamples << " " << fIfft[60]/numSamples;
//delete allocated resources
        return a.exec();
    }

Answer 2

In that case, I suspect it's a buffer-overrun. What does fileSize come back with in all cases?

Are you sure it's ios::end and not ios_base::end? Does it make a difference?

---

All you are doing in InverseFFT is taking the complex conjugate of the input data and then applying the FFT. That will not be enough to perform the inverse FFT.

See this page for details about how to perform the inverse FFT using just a forward FFT.