0
$\begingroup$

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

$\endgroup$

2 Answers 2

1
$\begingroup$

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();
    }
$\endgroup$
1
  • $\begingroup$ Good to know! Thanks for posting your answer. $\endgroup$
    – Peter K.
    Jul 14, 2017 at 13:39
0
$\begingroup$

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.

$\endgroup$
8
  • $\begingroup$ It works correctly for FFT/IFFT lengths of up to 2048 as I've said in my post. Also, I have just implemented the code in the dsp guide, which should presumably work. $\endgroup$
    – dsp_user
    Jul 10, 2017 at 13:14
  • $\begingroup$ Well, code doesn't implement the IFFT. Until it does, all bets are off. $\endgroup$
    – Peter K.
    Jul 10, 2017 at 13:24
  • 1
    $\begingroup$ I'll try some of the methods found in the link you provided and report back with the results. Thanks $\endgroup$
    – dsp_user
    Jul 10, 2017 at 13:41
  • $\begingroup$ Still no luck. I've updated my answer. $\endgroup$
    – dsp_user
    Jul 11, 2017 at 7:20
  • $\begingroup$ It makes no difference whether we use ios_base::end or ios::end. The input sound file (samples) are read correctly (I've used a hex editor and manually checked the values). fileSize is 32KB as expected. $\endgroup$
    – dsp_user
    Jul 11, 2017 at 11:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.