IFFT to Wavetable ?

Hi guys,

I'm trying to implement an algorithm where you fill in real and imaginary 
arrays then do an Inverse FFT to convert them to a wavetable.

The code I'm basing it from is like this:

double real[4];
double imag[4];

// fill them in

double wavetable[8];

IFFT(real,imag,wavetable);

It's using an FFT library (in a DLL) so I don't know what's going on. The 
thing I don't understand is how I get a real array of length 8 from a 
complex array of length 4. The FFT I use takes real and imag arrays and 
returns the result in the same arrays, like this:

FFT(int dir, int n,double* real,double* imag);

...so how do I combine them to make the wavetable ? The wavetable is 
supposed to be perfectly loop-able also.

Thanks.

"Robert A." <invalid@invalid.org> wrote in message 
news:WBT4g.168$UY5.152@fe10.lga...
> Hi guys,
>
> I'm trying to implement an algorithm where you fill in real and imaginary 
> arrays then do an Inverse FFT to convert them to a wavetable.
>
> The code I'm basing it from is like this:
>
> double real[4];
> double imag[4];
>
> // fill them in
>
> double wavetable[8];
>
> IFFT(real,imag,wavetable);
>
> It's using an FFT library (in a DLL) so I don't know what's going on. The 
> thing I don't understand is how I get a real array of length 8 from a 
> complex array of length 4.
Spect you are getting 4 real values and 4 imag values returned - you want to 
find how they are  ordered ?

> The FFT I use takes real and imag arrays and returns the result in the 
> same arrays, like this:
>
> FFT(int dir, int n,double* real,double* imag);
>
> ...so how do I combine them to make the wavetable ? The wavetable is 
> supposed to be perfectly loop-able also.
>
Why not use some other freely available fft routine to generate some test 
vectors and results then you can compare your outputs to find out what is 
real/imaginary and how ordered?

Best of luck - Mike

i'm gratified to see someone using the term "wavetable" as it was
originally used for "wavetable synthesis" and not in the sense of
"wavetable" soundcards or sampling keyboards.

Robert A. wrote:
> I'm trying to implement an algorithm where you fill in real and imaginary
> arrays then do an Inverse FFT to convert them to a wavetable.
>
> The code I'm basing it from is like this:
>
> double real[4];
> double imag[4];
>
> // fill them in
>
> double wavetable[8];
>
> IFFT(real,imag,wavetable);
>
> It's using an FFT library (in a DLL) so I don't know what's going on. The
> thing I don't understand is how I get a real array of length 8 from a
> complex array of length 4.

conjugate symmetry is implied. if you compute the DFT (which is what
the FFT is) of a purely real sequece of length N, you get this symmetry

     X[N-k] = conj{ X[k] }

which means X[0] and X[N/2] are purely real (have zero imaginary parts)
and the N/2-1 complex numbers in the 2nd half are complex conjugates of
those in the 1st have.  so there are only N unique real numbers coming
out (as there are N real numbers going in) of the FFT.

if there are only N/2 complex numbers used from X[0] to X[N/2 -1], then
usually what they do is put the real part of X[N/2] (the imaginary part
we *know* is zero) into the imaginary part of X[0] (which really should
be zero, but they're not wasting the space).

> The FFT I use takes real and imag arrays and
> returns the result in the same arrays, like this:
>
> FFT(int dir, int n,double* real,double* imag);
>
> ...so how do I combine them to make the wavetable ? The wavetable is
> supposed to be perfectly loop-able also.

you know what (harmonic) sinusoidal components you have, right?  do you
have them in terms of phase and magnitude?  you'll have to expand them
to real and imaginary parts (polar to rectangular conversion).  then
stick the real and imaginary parts into the real and imag arrays.

r b-j

Robert A. wrote:
> Hi guys,
>
> I'm trying to implement an algorithm where you fill in real and imaginary
> arrays then do an Inverse FFT to convert them to a wavetable.
>
> The code I'm basing it from is like this:
>
> double real[4];
> double imag[4];
>
> // fill them in
>
> double wavetable[8];
>
> IFFT(real,imag,wavetable);
>
> It's using an FFT library (in a DLL) so I don't know what's going on. The
> thing I don't understand is how I get a real array of length 8 from a
> complex array of length 4. The FFT I use takes real and imag arrays and
> returns the result in the same arrays, like this:
>
> FFT(int dir, int n,double* real,double* imag);
>
> ...so how do I combine them to make the wavetable ?

A standard IFFT, when used to produce a purely real result
vector (minus rounding error) requires redundant inputs.
So you need double real[8], imag[8], with all of the real
elements duplicated and all of the imaginary elements
duplicated with negation.

This input redundancy is built into standard IFFT routines
because they are designed for result vectors that are complex
as well as real.

Google "FFT Hermitian symmetry".


IMHO. YMMV.
-- 
rhn A.T nicholson d.0.t C-o-M

obert A. wrote:
> Hi guys,

> I'm trying to implement an algorithm where you fill in real and imaginary
> arrays then do an Inverse FFT to convert them to a wavetable.

> The code I'm basing it from is like this:

> double real[4];
> double imag[4];

> // fill them in

> double wavetable[8];

> IFFT(real,imag,wavetable);

> It's using an FFT library (in a DLL) so I don't know what's going on. The
> thing I don't understand is how I get a real array of length 8 from a
> complex array of length 4. The FFT I use takes real and imag arrays and
> returns the result in the same arrays, like this:

> FFT(int dir, int n,double* real,double* imag);

> ...so how do I combine them to make the wavetable ?

A standard IFFT, when used to produce a purely real result
vector (minus rounding error), requires redundant inputs.
So you need double real[8], imag[8]; with all of the real
elements (minus the DC and Nyquist terms) duplicated as in
a mirror image, and all of the imaginary elements duplicated
with negation,

This input redundancy is built into standard IFFT routines
because they are designed for result vectors which can be
complex as well as real, and thus require 2N degrees of
freedom.

Google "FFT Hermitian symmetry".

IMHO. YMMV.
-- 
rhn A.T nicholson d.0.t C-o-M

Hello,

After you explained it to me I think I understand now but I'm not positive 
yet.

I also wasn't converting from polar to rectangular.

So is the following what I need to do ?

***

double amplitude[4];
double phase[4];

// fill them in

// now make new arrays for the FFT (twice the size)

double real[8];
double imag[8];

// convert first half to rectangular and fill in part of the array

for(i=0;i<4;i++){
 real[i]=amplitude[i]*cos(phase[i]);
 imag[i]=amplitude[i]*sin(phase[i]);
}

// now set the second half to the complex conjugates

for(i=0;i<4;i++){
 real[i+4]=real[4-i];
 imag[i+4]=-imag[4-i];
}

// I think that part might be messed up because of element 0 and element 4 ?

// then do the inverse FFT

FFT(-1,8,real,imag);

// then my wavetable is the real components ?

double wavetable[8];

for(i=0;i<8;i++){
 wavetable[i]=real[i];
}

// probably have to scale the wavetable to the range [-1.0,1.0]

***

Am I on the right track here ? I think I almost have it. Thanks for helping 
me out.

Take care,
Robert A.



"robert bristow-johnson" <rbj@audioimagination.com> wrote in message 
news:1146419265.590185.155880@i39g2000cwa.googlegroups.com...
>
> i'm gratified to see someone using the term "wavetable" as it was
> originally used for "wavetable synthesis" and not in the sense of
> "wavetable" soundcards or sampling keyboards.
>
> Robert A. wrote:
>> I'm trying to implement an algorithm where you fill in real and imaginary
>> arrays then do an Inverse FFT to convert them to a wavetable.
>>
>> The code I'm basing it from is like this:
>>
>> double real[4];
>> double imag[4];
>>
>> // fill them in
>>
>> double wavetable[8];
>>
>> IFFT(real,imag,wavetable);
>>
>> It's using an FFT library (in a DLL) so I don't know what's going on. The
>> thing I don't understand is how I get a real array of length 8 from a
>> complex array of length 4.
>
> conjugate symmetry is implied. if you compute the DFT (which is what
> the FFT is) of a purely real sequece of length N, you get this symmetry
>
>     X[N-k] = conj{ X[k] }
>
> which means X[0] and X[N/2] are purely real (have zero imaginary parts)
> and the N/2-1 complex numbers in the 2nd half are complex conjugates of
> those in the 1st have.  so there are only N unique real numbers coming
> out (as there are N real numbers going in) of the FFT.
>
> if there are only N/2 complex numbers used from X[0] to X[N/2 -1], then
> usually what they do is put the real part of X[N/2] (the imaginary part
> we *know* is zero) into the imaginary part of X[0] (which really should
> be zero, but they're not wasting the space).
>
>> The FFT I use takes real and imag arrays and
>> returns the result in the same arrays, like this:
>>
>> FFT(int dir, int n,double* real,double* imag);
>>
>> ...so how do I combine them to make the wavetable ? The wavetable is
>> supposed to be perfectly loop-able also.
>
> you know what (harmonic) sinusoidal components you have, right?  do you
> have them in terms of phase and magnitude?  you'll have to expand them
> to real and imaginary parts (polar to rectangular conversion).  then
> stick the real and imaginary parts into the real and imag arrays.
>
> r b-j
> 

Robert A. wrote:
> Hello,
>
> After you explained it to me I think I understand now but I'm not positive
> yet.

you're getting very close.  here is how i would do it:


// size of wavetable.  needs to be even, maybe a power or 2.
#define N 128



double amplitude[N/2];
double phase[N/2];


// fill them in...

for (int k=1; k<N/2; k++)
    {
    amplitude[k] = ...;     // amplitude of kth harmonic
    phase[k] = ...;    // phase of kth harmonic
    }

phase[0] = 0.0;       // this could be +/- pi or 0, but nothing else

 // my suggestion is that amplitude[N/2] be zero

// now make new arrays for the FFT (twice the size)
double real[N];
double imag[N];

// convert first half to rectangular and fill in part of the array


double scaler = (double)N/2.0;
 // maybe just scaler = 1.0/2.0 if iFFT doesn't scale by 1/N

real[0] = 2.0*scaler*amplitude[k];   // DC and Nyquist are a little
different
imag[0] = 0.0;     // just to make awul damn sure.

double energy = rea[0]*real[0];

for (int k=1; k<N/2; k++)
    {
    real[k] = scaler*amplitude[k]*cos(phase[k]);
    imag[k] = scaler*amplitude[k]*sin(phase[k]);

    energy += rea[k]*real[k] + imag[k]*imag[k];
    }

// now set the second half to the complex conjugates

real[N/2] = 0.0;   // amplitude[N/2] = zero
imag[N/2] = 0.0;    // just to make awul damn sure.

for (int k=N/2+1; k<N; k++)
    {
    real[k] = real[N-k];
    imag[k] = -imag[N-k];
    }


// I think that part might be messed up
// because of element 0 and element N/2 ?

    //  it was, but i fixed it.

// then do the inverse FFT

FFT(-1, N, real, imag);
   // i am presuming that iFFT has a 1/N in it.

// then my wavetable is the real components ?
    //   yup. the imag components should all be zero.

double wavetable[N];
double imag_error_limit = log((double)N)*energy*1e-15;
int imag_error = 0;

for (int n=0; n<N; n++)
    {
    wavetable[n] = real[n];
    if (imag[n]*imag[n] > imag_error_limit)
        {
        imag_error++;
        }
    }

if (imag_error)
    {
    printf("i did something wrong.  imag_error = %d", imag_error);
    }

// data is in wavetable[] and imag_error should be 0.
// scaling of wavetable[] can be done properly by scaling amplitude[]
correctly.


so, you had the right idea. does this do it for you?

BTW, i forgot to mention that the DFT (which the FFT is a fast
implementation of) maps a discrete-time periodic function of period N
to another discrete-frequency periodic function of period N.  so the
result is "perfectly loop-able".  if you want to get more anal
wavetable mathematics than you would ever want, take a look at:

    http://www.musicdsp.org/files/Wavetable-101.pdf

r b-j

Hi Robert,

It's definitely working great. I just kind of copied and pasted the code for 
now but I intend to really go over it over so I can learn exactly what's 
going on. I really appreciate you taking the time to help me out.

One thing I noticed when working on this stuff is that the FFT routine in 
Numerical Recipes in C gives opposite signs for the imaginary numbers, I 
always considered that book the standard and I'm not sure what to think now. 
I'm comparing the results to the example in Understand Digital Signal 
Processing by Lyons (I just got that book)

I switched over to this routine (the FFT one at the bottom)

http://astronomy.swin.edu.au/~pbourke/other/dft/

...it's a lot easier to work with and the numbers match the Lyons book.

Thanks again.

PS - I'm sorry it took so long to respond


"robert bristow-johnson" <rbj@audioimagination.com> wrote in message 
news:1146625438.510983.281330@e56g2000cwe.googlegroups.com...
> Robert A. wrote:
>> Hello,
>>
>> After you explained it to me I think I understand now but I'm not 
>> positive
>> yet.
>
> you're getting very close.  here is how i would do it:
>
>
> // size of wavetable.  needs to be even, maybe a power or 2.
> #define N 128
>
>
>
> double amplitude[N/2];
> double phase[N/2];
>
>
> // fill them in...
>
> for (int k=1; k<N/2; k++)
>    {
>    amplitude[k] = ...;     // amplitude of kth harmonic
>    phase[k] = ...;    // phase of kth harmonic
>    }
>
> phase[0] = 0.0;       // this could be +/- pi or 0, but nothing else
>
> // my suggestion is that amplitude[N/2] be zero
>
> // now make new arrays for the FFT (twice the size)
> double real[N];
> double imag[N];
>
> // convert first half to rectangular and fill in part of the array
>
>
> double scaler = (double)N/2.0;
> // maybe just scaler = 1.0/2.0 if iFFT doesn't scale by 1/N
>
> real[0] = 2.0*scaler*amplitude[k];   // DC and Nyquist are a little
> different
> imag[0] = 0.0;     // just to make awul damn sure.
>
> double energy = rea[0]*real[0];
>
> for (int k=1; k<N/2; k++)
>    {
>    real[k] = scaler*amplitude[k]*cos(phase[k]);
>    imag[k] = scaler*amplitude[k]*sin(phase[k]);
>
>    energy += rea[k]*real[k] + imag[k]*imag[k];
>    }
>
> // now set the second half to the complex conjugates
>
> real[N/2] = 0.0;   // amplitude[N/2] = zero
> imag[N/2] = 0.0;    // just to make awul damn sure.
>
> for (int k=N/2+1; k<N; k++)
>    {
>    real[k] = real[N-k];
>    imag[k] = -imag[N-k];
>    }
>
>
> // I think that part might be messed up
> // because of element 0 and element N/2 ?
>
>    //  it was, but i fixed it.
>
> // then do the inverse FFT
>
> FFT(-1, N, real, imag);
>   // i am presuming that iFFT has a 1/N in it.
>
> // then my wavetable is the real components ?
>    //   yup. the imag components should all be zero.
>
> double wavetable[N];
> double imag_error_limit = log((double)N)*energy*1e-15;
> int imag_error = 0;
>
> for (int n=0; n<N; n++)
>    {
>    wavetable[n] = real[n];
>    if (imag[n]*imag[n] > imag_error_limit)
>        {
>        imag_error++;
>        }
>    }
>
> if (imag_error)
>    {
>    printf("i did something wrong.  imag_error = %d", imag_error);
>    }
>
> // data is in wavetable[] and imag_error should be 0.
> // scaling of wavetable[] can be done properly by scaling amplitude[]
> correctly.
>
>
> so, you had the right idea. does this do it for you?
>
> BTW, i forgot to mention that the DFT (which the FFT is a fast
> implementation of) maps a discrete-time periodic function of period N
> to another discrete-frequency periodic function of period N.  so the
> result is "perfectly loop-able".  if you want to get more anal
> wavetable mathematics than you would ever want, take a look at:
>
>    http://www.musicdsp.org/files/Wavetable-101.pdf
>
> r b-j
> 

Robert A. wrote:

...
> One thing I noticed when working on this stuff is that the FFT routine in
> Numerical Recipes in C gives opposite signs for the imaginary numbers, I
> always considered that book the standard and I'm not sure what to think now.
> I'm comparing the results to the example in Understand Digital Signal
> Processing by Lyons (I just got that book)

There is some disagreement in the sciences to what actually the Fourier
transform is. For the continuous case look at the discussion here,

http://mathworld.wolfram.com/FourierTransform.html

following equations (9) and (10). Some of this confusion carries over
into the discrete world, where physicists and engineers have conjugate
definitions of the transfrom (Numerical Recipes takes the physicists
stand, pfui!). On top of that, sometimes it is useful to have different
scaling factor in front of the transform. This often causes confusion
when you need to go there and back, and want to land in the starting
position (for example in frequency domain filtering). The archives of
comp.dsp are awash with sorry souls that have shipwrecked in the
treacherous waters of the DFT.

Regards,
Andor

"Andor" <andor.bariska@gmail.com> wrote in message 
news:1147169896.111443.32510@y43g2000cwc.googlegroups.com...
> Robert A. wrote:
>
> ...
>> One thing I noticed when working on this stuff is that the FFT routine in
>> Numerical Recipes in C gives opposite signs for the imaginary numbers, I
>> always considered that book the standard and I'm not sure what to think 
>> now.
>> I'm comparing the results to the example in Understand Digital Signal
>> Processing by Lyons (I just got that book)
>
> There is some disagreement in the sciences to what actually the Fourier
> transform is. For the continuous case look at the discussion here,
>
> http://mathworld.wolfram.com/FourierTransform.html
>
> following equations (9) and (10). Some of this confusion carries over
> into the discrete world, where physicists and engineers have conjugate
> definitions of the transfrom (Numerical Recipes takes the physicists
> stand, pfui!).


I'm sure the physicists use the version that Fourier came up with to solve 
physics problems. He first applied it to heat flow.

Clay