FFTw - Manually regenerating the original signal

Hi!

I'm totally new to FT:s (obviosly) and don't know much of the math involved.

I have a real input signal (Fs=44100) which I've transformed to a 
complex and I was hoping to be able to regenerate the original signal 
manually - just to get the hang of it. It comes "close" :-)


I'm using these fftwl_plans:

long double *in;
fftwl_complex *out;
fftwl_plan plan_forward, plan_backward;

in = (long double*) fftwl_malloc(sizeof(long double)*fft_length);
out = fftwl_alloc_complex(fft_length);

plan_forward = fftwl_plan_dft_r2c_1d( fft_length, in, out, FFTW_EXHAUSTIVE);

plan_backward = fftwl_plan_dft_c2r_1d( fft_length, out, in, 
FFTW_EXHAUSTIVE);


I can transform forward and then backward and get the original signal 
back (after normalizing it with 1/fft_length).

fft_length=10 bins=fft_length/2+1=6

The forward transformation gives me this (unnormalized):

        binFreq ->       real       img  ->        amp     phase
       0.000000 -> {-0.398685, 0.000000} -> { 0.398685, 3.141593}
    4410.000000 -> {-1.017359, 0.591153} -> { 1.176640, 2.615212}
    8820.000000 -> { 1.187429,-0.732900} -> { 1.395396,-0.552982}
   13230.000000 -> { 0.515134,-1.726352} -> { 1.801570,-1.280813}
   17640.000000 -> { 0.406585,-0.089637} -> { 0.416349,-0.216991}
   22050.000000 -> { 0.383167, 0.000000} -> { 0.383167, 0.000000}

Here's the input signal, the signal after forward+backward FT 
(+normalization in the time domain) and then my two generated signals - 
one generated from the complex numbers and one from amplitude and phase. 
The error in my generated signals seems to be the same.

smp         in   fwd&back  gen(cplx) gen(AmpPhase)
   0   0.216806   0.216806   0.215254   0.215254
   1   0.141786   0.141786   0.063601   0.063601
   2  -0.561061  -0.561061  -0.562613  -0.562613
   3  -0.483457  -0.483457  -0.561642  -0.561642
   4   0.311432   0.311432   0.309880   0.309880
   5   0.341063   0.341063   0.262877   0.262877
   6   0.093564   0.093564   0.092012   0.092012
   7   0.285536   0.285536   0.207351   0.207351
   8  -0.068499  -0.068499  -0.070051  -0.070051
   9  -0.675855  -0.675855  -0.754040  -0.754040


Here's some C++ code describing how I try to do it:

typedef struct {
   long double amp;
   long double phase;
} AmpPhase;

inline void complex2AmpPhase( fftwl_complex &cpl, AmpPhase &ap ) {
   long double r,i;
   r=cpl[0];
   i=cpl[1];
   ap.amp = sqrtl( r*r + i*i );
   ap.phase = atan2l(i,r);
}

// binFunc( long bin, long sample ) =
//  binFreq(bin) * sample * M_TWO_PI / Fs
// ...but is somewhat optimized here

inline long double binFunc( long bin, long sample ) {
     return bin * sample * M_TWO_PI/fft_length;
}

// regenerate the signal from the unnormalized FT

for(int samp=0;samp<fft_length;samp++) {
   gen_cpx[samp]=0.0;
   gen_ap[samp]=0.0;
}

for( int b=0; b<bins; b++ ) {
   complex2AmpPhase( out[b], ap );

   for(int samp=0;samp<fft_length;samp++) {

     gen_cpx[samp] +=
       out[b][0] * cosl( binFunc(b,samp) ) -
       out[b][1] * sinl( binFunc(b,samp) );

     gen_ap[samp] += ap.amp * cosl( binFunc(b,samp) + ap.phase );
   }
}


As I said, I normalize the backwards transformation with 1/fft_length 
and get the correct signal back but I need to normalize my generated 
signals with 2/fft_length to get as close as I'm getting so I must have 
messed up bigtime somewhere. :)

If I instead normalize in the frequency domain, the backwards FT still 
gives me the correct answer but then my generated signals gets half the 
amplitude (same normalization problem).

If I let my insignal be a clean sin() spot on one of the bin frequencies 
(8820), this is what I get after normalization:

        binFreq ->       real       img  ->        amp     phase
       0.000000 -> { 0.000000, 0.000000} -> { 0.000000, 0.000000}
    4410.000000 -> { 0.000000, 0.000000} -> { 0.000000, 0.865899}
    8820.000000 -> {-0.000000,-0.500000} -> { 0.500000,-1.570796}
   13230.000000 -> { 0.000000,-0.000000} -> { 0.000000,-0.994491}
   17640.000000 -> { 0.000000,-0.000000} -> { 0.000000,-0.484800}
   22050.000000 -> { 0.000000, 0.000000} -> { 0.000000, 0.000000}

And if I let it be a clean sin() between two bins (11025):

        binFreq ->       real       img  ->        amp     phase
       0.000000 -> { 0.100000, 0.000000} -> { 0.100000, 0.000000}
    4410.000000 -> { 0.123607,-0.000000} -> { 0.123607,-0.000000}
    8820.000000 -> { 0.323607,-0.000000} -> { 0.323607,-0.000000}
   13230.000000 -> {-0.323607, 0.000000} -> { 0.323607, 3.141593}
   17640.000000 -> {-0.123607, 0.000000} -> { 0.123607, 3.141593}
   22050.000000 -> {-0.100000, 0.000000} -> { 0.100000, 3.141593}

And in the generated signals I get from the above, it's clear that 
something fishy is going on:

smp         in   fwd&back  gen(cplx) gen(AmpPhase)
   0   0.000000  -0.000000  -0.000000  -0.000000
   1   1.000000   1.000000   0.600000   0.600000
   2   0.000000   0.000000   0.000000   0.000000
   3  -1.000000  -1.000000  -0.400000  -0.400000
   4  -0.000000  -0.000000  -0.000000  -0.000000
   5   1.000000   1.000000   0.600000   0.600000
   6   0.000000   0.000000   0.000000   0.000000
   7  -1.000000  -1.000000  -0.400000  -0.400000
   8  -0.000000  -0.000000  -0.000000  -0.000000
   9   1.000000   1.000000   0.600000   0.600000

Sorry for the long post. Any help would be much appreciated!

Kind regards,
Ted Lyngmo

Ted Lyngmo wrote:
> Hi!
<snip>
>
> If I instead normalize in the frequency domain, the backwards FT still
> gives me the correct answer but then my generated signals gets half the
> amplitude (same normalization problem).
>

I think that this is as good as it gets. May be worth a trawl
through the documentation of FFTw concerning normalization - it's
sort of *your* problem as a coder using FFTw as I recall :)

I do remember being able to reproduce the original PCM signal by
round-tripping through FFTw - fftw_plan_dft_1d(len, h, H, FFTW_FORWARD,
FFTW_ESTIMATE) followed by fftw_plan_dft_1d(length, H, h, FFTW_BACKWARD, 
FFTW_ESTIMATE).

<snip>
> Kind regards,
> Ted Lyngmo

-- 
Les Cargill

Hi and thanks for your answer Les!

Les Cargill wrote:
> Ted Lyngmo wrote:
>>
>> If I instead normalize in the frequency domain, the backwards FT still
>> gives me the correct answer but then my generated signals gets half the
>> amplitude (same normalization problem).
>
> I think that this is as good as it gets. May be worth a trawl
> through the documentation of FFTw concerning normalization - it's
> sort of *your* problem as a coder using FFTw as I recall :)

4.8.2 The 1d Real-data DFT: "Like FFTW's complex DFT, these transforms 
are unnormalized. In other words, applying the real-to-complex (forward) 
and then the complex-to-real (backward) transform will multiply the 
input by n."

But I found this at the bottom of "4.8.3 1d Real-even DFTs" (which I'm 
not using):

"These definitions correspond directly to the unnormalized DFTs used 
elsewhere in FFTW (hence the factors of 2 in front of the summations)."

Perhaps I mixing stuff up when trying to reconstruct the original signal 
transformed with "The 1d Real-data DFT" by using cos & sin like this?

for(int b=0; b<bins; b++) {
   for(int samp=0;samp<fft_length;samp++) {
     gen_cpx[samp] +=
       out[b][0] * cosl( b*samp*M_TWO_PI/fft_length ) -
       out[b][1] * sinl( b*samp*M_TWO_PI/fft_length );
   }
}

> I do remember being able to reproduce the original PCM signal by
> round-tripping through FFTw - fftw_plan_dft_1d(len, h, H, FFTW_FORWARD,
> FFTW_ESTIMATE) followed by fftw_plan_dft_1d(length, H, h, FFTW_BACKWARD,
> FFTW_ESTIMATE).

Yes, that works for me too (that's the "fwd&back" column in my examples).

in -> FT_forward -> FT_backward -> normalize (1/n) -> out
in = out

in -> FT_forward -> normalize (1/n) -> FT_backward -> out
in = out

Br,
Ted

I just noticed that my original message isn't displayed correctly when 
viewed online at http://www.dsprelated.com/showthread/comp.dsp/287119-1.php

The whole part where I showed how I regenerate the signal isn't 
displayed. :-/

It looks ok at google groups 
https://groups.google.com/forum/#!original/comp.dsp/E5Djuw9TwQ8/es_m7Wyk6oQJ

Br,
Ted

 > 
https://groups.google.com/forum/#!original/comp.dsp/E5Djuw9TwQ8/es_m7Wyk6oQJ 


I think I stumbled across how to regenerate the original signal by just 
looking at the numbers for long enough, but still don't get the math 
behind it. I've just compensated for what I observed.

This seems to work:

I need to adjust the regenerated signal with what's in the 0 Hz bin's 
real part.

And if I'm using an even fft length, I need to also adjust with the 
highest freqency bin's real part. Adding for odd samples and subtracting 
for even ones.


I'm trying to post code again even though dsprelated seems to mess it 
up. This regenerates the signal from the normalized FT that is in 'out':

for( int b=0; b < bins; b++ ) {
   for( int samp=0; samp < fft_length; samp++ ) {

     gen_cpx[samp] +=
       out[b][0] * cosl( binFunc(b,samp) ) -
       out[b][1] * sinl( binFunc(b,samp) );

   }
}

// normalize and adjust

for( int samp=0; samp<gen_length; samp++ ) {
   gen_cpx[samp] *= 2.0;
   gen_cpx[samp] -= out[0][0]; // offset

   if( !(fft_length&1) ) { // only even fft lenghts
     if( samp&1 ) { // odd sample
       gen_cpx[samp] += out[bins-1][0];
     } else {
       gen_cpx[samp] -= out[bins-1][0];
     }
   }
}


I've tried with different fft_lengths and signals (and noise) and so far 
it seems to regenerate the signal perfectly.

Br,
Ted

 > 
https://groups.google.com/forum/#!original/comp.dsp/E5Djuw9TwQ8/es_m7Wyk6oQJ

I think I stumbled across how to regenerate the original signal by just 
looking at the numbers for long enough, but still don't get the math 
behind it. I've just compensated for what I observed.

This seems to work:

I need to adjust the regenerated signal with what's in the 0 Hz bin's 
real part.

And if I'm using an even fft length, I need to also adjust with the 
highest freqency bin's real part. Adding for odd samples and subtracting 
for even ones.


I'm trying to post code again even though dsprelated seems to mess it 
up. This regenerates the signal from the normalized FT that is in 'out':

for( int b=0; b < bins; b++ ) {
   for( int samp=0; samp < fft_length; samp++ ) {

     gen_cpx[samp] +=
       out[b][0] * cosl( binFunc(b,samp) ) -
       out[b][1] * sinl( binFunc(b,samp) );

   }
}

// normalize and adjust

for( int samp=0; samp < gen_length; samp++ ) {
   gen_cpx[samp] *= 2.0;
   gen_cpx[samp] -= out[0][0]; // offset

   if( !(fft_length&1) ) { // only even fft lenghts
     if( samp&1 ) { // odd sample
       gen_cpx[samp] += out[bins-1][0];
     } else {
       gen_cpx[samp] -= out[bins-1][0];
     }
   }
}


I've tried with different fft_lengths and signals (and noise) and so far 
it seems to regenerate the signal perfectly.

Br,
Ted

>I just noticed that my original message isn't displayed correctly when 
>viewed online at
http://www.dsprelated.com/showthread/comp.dsp/287119-1.php
>
>The whole part where I showed how I regenerate the signal isn't 
>displayed. :-/

Indeed!  Sorry about this. I think I found and corrected the error in the
code and it should look better now. Never hesitate to report problems like
this to me through the contact form on DSPRelated.com.  I always
appreciate any feedback that help me improve the website.  Thanks.

---------------------------------------
Posted through http://www.DSPRelated.com