swapping realft() for fftw, need help

Hi,

I have been refactoring some old SMS code that uses a non-GPL'ed (and
probably much slower) fft function from a book called "Numerical Recipes in
C." One of the first things I am doing is changing this to use fftw3f... I
started by focusing on the resynthesis, which uses the inverse-fft method
for synthesizing a bank of sinewaves. 

But I have run into a brick wall: I thought all would have to do is pack
an array of

(float) [real0,imag0, real1, imag1, ..., realN, imagN]


into a complex array of:
[ [0][real,imag], [1][real,imag], ... [N][real,imag] ]  //sorry, this is
bad notation, but hopefully you get the point

but when I run the c2r fft on this, I get junk. Well, it isn't junk,
because there are periodic distortions, but I can't figure out what they
mean.  There is also a big magnitude problem going on, my waveform's
magnitude is 300x that of the realft transform! Lastly, there is a strange
hole in my waveform... or two actually, and they are consistent.  So,
hopefully someone with experience can see what I am doing wrong, or suggest
ways to further debug this.

Here are some plots of the output waveforms from each of the inverse fft
algorithms.
- horribly distorted fftw on top, realft on bottom:
http://www.teafordragons.com/rte/img/sms/fftComparison.png
- same as above, but only 256 samples:
http://www.teafordragons.com/rte/img/sms/fftComparison-zoomed.png
- just the messed up fftw waveform, showing the periodic 'gaps':
http://www.teafordragons.com/rte/img/sms/fftw-waveform-zoomed.png
- Hide quoted text -


The version of "Numerical Recipes in C" used is now obsolete and is online
at http://www.nrbook.com/a/bookcpdf.php
Chapters 12.2 and 12.3 contain the important code that I am trying to
replace, which to me looks identical to the fftw explanations, sadly.

The following is a summary of the important parts of the code and I can
supply anything else that is necessary:

------------------------------------------------------------------
/* make fftw plan and buffer for Inverse real FFT (complex->real)*/
int initInverseFFTW( SYNTH_PARAMS *pSynthParams)
{
        int sizeFFT = pSynthParams->sizeHop ;  
        pSynthParams->pCfftIn =  fftwf_malloc(sizeof(fftwf_complex) *
(sizeFFT / 2 + 1));
        pSynthParams->pFfftOut = fftwf_malloc(sizeof(float) * sizeFFT);
        pSynthParams->fftPlan =  fftwf_plan_dft_c2r_1d( sizeFFT,
pSynthParams->pCfftIn,
                                                     
pSynthParams->pFfftOut, FFTW_ESTIMATE);
        
        /*debugging realft */
        pSynthParams->realftOut = (float *) calloc((sizeFFT<<1)+1,
sizeof(float));
        return(1);
} 


/* synthesis of one frame of the deterministic component using the IFFT
*/
static int SineSynthIFFT (SMS_DATA *pSmsData, float *pFBuffer,
                          SYNTH_PARAMS *pSynthParams)
{
        long sizeFft = pSynthParams->sizeHop << 1;
        long iHalfSamplingRate = pSynthParams->iSamplingRate >> 1;
        long sizeMag = pSynthParams->sizeHop;
        long nBins = 8;
        long  nTraj = pSmsData->nTraj;
        long iFirstBin, k, i, l, b;
        float fMag=0, fFreq=0, fPhase=0, fLoc, fSin, fCos, fBinRemainder,
                fTmp, fNewMag,  fIndex;

//        pFFftBuffer = (float *) calloc(sizeFft+1, sizeof(float));
        memset (pSynthParams->realftOut, 0, (sizeFft +1) *
sizeof(float));
        for (i = 0; i < nTraj; i++)
        {
                if (((fMag = pSmsData->pFMagTraj[i]) > 0) &&
                    ((fFreq = pSmsData->pFFreqTraj[i]) <
iHalfSamplingRate))
                {
                        if (pSynthParams->previousFrame.pFMagTraj[i] <=
0)
                                pSynthParams->previousFrame.pFPhaTraj[i]
=
                                        TWO_PI * ((random() - HALF_MAX) /
HALF_MAX);
                        // can fMag here be stored as magnitude instead of
DB within smsData?
                        fMag = TO_MAG (fMag);
                        fTmp = pSynthParams->previousFrame.pFPhaTraj[i] +
                                (TWO_PI * fFreq /
pSynthParams->iSamplingRate) * sizeMag;
                        fPhase = fTmp - floor(fTmp / TWO_PI) * TWO_PI;
                        fLoc = sizeFft * fFreq  /
pSynthParams->iSamplingRate;
                        iFirstBin = (int) fLoc - 3;
                        fBinRemainder = fLoc - floor (fLoc);
                        fSin = SinTab (fPhase);
                        fCos = SinTab (fPhase + PI_2);
                        for (k = 1, l = iFirstBin; k <= nBins; k++, l++)
                        {
                                fIndex = (k - fBinRemainder);
                                if (fIndex > 7.999) fIndex = 0;
                                fNewMag = fMag * SincTab (fIndex);
                                if (l > 0 && l < sizeMag)
                                {
                                        pSynthParams->realftOut[l*2+1] +=
fNewMag * fCos;
                                        pSynthParams->realftOut[l*2] +=
fNewMag * fSin;
                                }
                                else if (l == 0)
                                {
                                        pSynthParams->realftOut[0] += 2 *
fNewMag * fSin;
                                }
                                else if (l < 0)
                                {
                                        b = abs(l);
                                        pSynthParams->realftOut[b*2+1] -=
fNewMag * fCos;
                                        pSynthParams->realftOut[b*2] +=
fNewMag * fSin;
                                }
                                else if (l > sizeMag)
                                {
                                        b = sizeMag - (l - sizeMag);
                                        pSynthParams->realftOut[b*2+1] -=
fNewMag * fCos;
                                        pSynthParams->realftOut[b*2] +=
fNewMag * fSin;
                                }
                                else if (l == sizeMag)
                                {
                                        pSynthParams->realftOut[1] += 2 *
fNewMag * fSin;
                                }
                        }
                }
                pSynthParams->previousFrame.pFMagTraj[i] = fMag;
                pSynthParams->previousFrame.pFPhaTraj[i] = fPhase;
                pSynthParams->previousFrame.pFFreqTraj[i] = fFreq;
        }

        // fill complex fftw buffer:
        // pSynthParams->realftOut[i*2] is complex (sin), while
fftw_complex[i][0] is real.
        // so they have to be switched -- I think this is backwards,
actually... ARG
        // shouldn't change anything that drastically in amplitude though
        for(i = 0; i < (sizeMag/2 +1); i++)
        {
                pSynthParams->pCfftIn[i][0] =
pSynthParams->realftOut[i*2];
                pSynthParams->pCfftIn[i][1] =
pSynthParams->realftOut[i*2+1];
        }
       // both versions of fft are done for now for testing
       fftwf_execute(pSynthParams->fftPlan);
       realft (pSynthParams->realftOut-1, sizeMag, -1);

#ifdef FFTW//not right yet...
        {
                //realft needs to be multiplied by 2/n, which is
apparently done elsewhere.
                //... so does this need to be mulitplied by 0.5 because
fftw has gain of 1/n ?
                // well since it is about 300x too big in amp, it isn't
the basis of the problem
                for(i = 0, k = sizeMag; i < sizeMag; i++, k++)
                        pFBuffer[i] += pSynthParams->pFfftOut[k] *
pSynthParams->pFDetWindow[i];
                for(i= sizeMag, k = 0; i < sizeFft; i++, k++)
                        pFBuffer[i] +=  pSynthParams->pFfftOut[k] *
pSynthParams->pFDetWindow[i];
        }
#else
        {
           
               
                for(i = 0, k = sizeMag; i < sizeMag; i++, k++)
                        pFBuffer[i] += pSynthParams->realftOut[k] *
pSynthParams->pFDetWindow[i];
                for(i= sizeMag, k = 0; i < sizeFft; i++, k++)
                        pFBuffer[i] +=  pSynthParams->realftOut[k] *
pSynthParams->pFDetWindow[i];

        }
#endif


        return (1);
}

------------------------------------------------------------------


This is my first time using FFTW3, so I wasn't expecting everything to go
smoothly.  But I am a little stuck... thanks to anyone who can help!

regards,
rich