IFFT to Wavetable ?

I found a small problem with the code:

"real[0] = 2.0*scaler*amplitude[k];   // DC and Nyquist are a little 
different"

what is k in that ? I think it's (n/2-1) but I'm not really sure.

Thanks.





"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:
> I found a small problem with the code:
>
> "real[0] = 2.0*scaler*amplitude[k];   // DC and Nyquist are a little different"
>
> what is k in that ? I think it's (n/2-1) but I'm not really sure.

oooops.  it's 0.  real[0] is the real part of the DC componet and
amplitude[0] is the amplitude of the 0th harmonic, or the DC component.
 dunno what "n" is.

it was from copy/paste and modify and i didn't do it perfectly right.

for musical waveforms, if i were you i would set the DC and Nyquist
components to zero and then there is no confusion with them.
particularly Nyquist (N/2) should be set explicitly to zero. save
yourself headache.

r b-j