FFTW3 fftw_plan_dft_r2c_1d and its inverse question

Hello,

I am running FreeBSD. I installed FFTW3 from the ports. I am not a native
english speaker so I apolozige for the mistakes.

The goal of my program is to perform the FFT in order to extract all the
required informations from the audio signal (frequencies, amplitudes,
phases) , same them in a file , and then rebuild the original audio signal
from the file.

A)
I try the FFTW c2r to go "backwards" (full source code at the end)

 plan_forward = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);
 plan_backward = fftw_plan_dft_c2r_1d(N, out, in, FFTW_ESTIMATE |
FFTW_PRESERVE_INPUT);

 fftw_execute(plan_forward);
 printf("\tDONE! (1= %f, 2= %f, 2000= %f, 20000= %f)\n", in[1], in[2],
in[2000], in[20000]);

 fftw_execute(plan_backward);
 printf("\tDONE! (1= %f, 2= %f, 2000= %f, 20000= %f)\n", in[1], in[2],
in[2000], in[20000]);

I don't get twice the same output.

That means that c2r didn't recreate the original signal from the result of
r2c. How to get the original one from the complex array?

B)
I can find the frequency and the amplitude
 frequency = i
 amplitude = (unsigned int)((2.0*sqrt(real*real + imag*imag)) / N);

but how can I find the phase?
(http://en.wikipedia.org/wiki/Frequency_spectrum say that it is possible
but not how)

C)
If I have a silence period of 0.2 seconds before starting a sinus of 120
Hz with an amplitude of 5000, will the silence be erased from the rebuilt
file?

Actually, in order to rebuild the original signal, I need to know the
frequency, the amplitude, the start time and the end time of each sinus?

D)
In my code, I had to put:
 if (ampli) {
  // fr�quence = i , amplitude = ampli
  ampli++;

to find the same amplitude as the one I put in gen_sound().
Why don't FFTW3 do the +1 itself?

---

Thanks for your help!

---

Here is the full source code:

#include <stdio.h>
#include <unistd.h>
#include "/usr/local/include/fftw3.h"	// THIS MAY BE WRONG ON YOUR
COMPUTER
#include <math.h>

#define SOUND_FILE	"nemo-gen.raw"	// Fichier entr&#4294967295;e
#define SOUND_OFILE	"nemo-gen.raw"	// Fichier sortie
#define CHANNELS	1		// Canaux (2=st&#4294967295;r&#4294967295;o, 1=mono)
#define SAMPLE_BYTES	2		// Nombre d'octets par sample (2=16bits, 1=8bits)
#define N		44100		// Fr&#4294967295;quence d'&#4294967295;chantillonage (pair)

#define N_OUT		N<<1 + 1
#define T_DELTA		1.0/N

// Generates 1 second of sound
void gen_sound() {
	double t=0;
	unsigned int i=0;
	short x;
	FILE *fp;

	fp=fopen(SOUND_OFILE, "w");
	if (!fp) return;

	while (i<N) {
		x=0;
		x += (short)(9000 * sin(2*M_PI * 122 * t));
		//x += (short)(9000 * sin(2*M_PI * 122 * t + M_PI)); // destructive
interference

		fwrite(&x, sizeof(short), 1, fp);

		t+=T_DELTA;
		i++;
	}

	fclose(fp);
}

int main(int argc, char* argv[]) {
	double *in;
	fftw_complex *out;
	fftw_plan plan_forward, plan_backward;
	FILE *fp;
	short sample[CHANNELS];
	unsigned int i, ampli;
	double real, imag;

	// G&#4294967295;n&#4294967295;ration du signal
	printf("Signal generation..\n");
	gen_sound();

	// Ouverture du fichier
	printf("Opening file..\n");
	fp=fopen(SOUND_FILE, "r");
	if (!fp) {
		printf("[ERREUR] Impossible d'ouvrir le fichier %s!\n", SOUND_FILE);
		return 1;
	}

	// Allocation de la m&#4294967295;moire
	printf("Allocating memory..\n");
	in = fftw_malloc(sizeof(double) * N);
	out = fftw_malloc(sizeof(fftw_complex) * N_OUT);

	// Cr&#4294967295;ation du plan pour le DFT
	printf("Creating plans..\n");
	plan_forward = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);
	plan_backward = fftw_plan_dft_c2r_1d(N, out, in, FFTW_ESTIMATE |
FFTW_PRESERVE_INPUT);

	// Chargement de l'&#4294967295;chantillon
	printf("Loading audio..\n");
	i=0;
	while (i<N) {
		fread(sample, sizeof(*sample), CHANNELS, fp);
		in[i]=sample[0]; // Left channel
		i++;
	}

	// Fermeture du fichier
	printf("Closing file..\n");
	fclose(fp);

	// Calcul de la transform&#4294967295;e
	printf("Executing DTF forward..\n");
	fftw_execute(plan_forward);
	printf("\tDONE! (1= %f, 2= %f, 2000= %f, 20000= %f)\n", in[1], in[2],
in[2000], in[20000]);

	// Calcul de la tranform&#4294967295;e inverse
	printf("Executing DTF backward..\n");
	fftw_execute(plan_backward);
	printf("\tDONE! (1= %f, 2= %f, 2000= %f, 20000= %f)\n", in[1], in[2],
in[2000], in[20000]);

	// On r&#4294967295;cup&#4294967295;re les harmoniques
	printf("Finding frequencies..\n");
	i=0;
	while (i<N_OUT) {
		real = out[i][0];
		imag = out[i][1];

		ampli = (unsigned int)((2.0*sqrt(real*real + imag*imag)) / N);

		if (ampli) {
			// fr&#4294967295;quence = i , amplitude = ampli
			ampli++;
			printf("Freq= %u, ampli= %u\n", i, ampli);
		}

		i++;
	}

	// On lib&#4294967295;re les ressources
	printf("Destroying resources..\n");
	fftw_destroy_plan(plan_forward);
	fftw_destroy_plan(plan_backward);
	fftw_free(in);
	fftw_free(out);

	return 0;
}



		
This message was sent using the Comp.DSP web interface on
www.DSPRelated.com

dinsoft wrote:
> That means that c2r didn't recreate the original signal from the result of
> r2c. How to get the original one from the complex array?

You need to divide by N; see question 3.9 in the FFTW FAQ.

>
> B)
> I can find the frequency and the amplitude
>  frequency = i
>  amplitude = (unsigned int)((2.0*sqrt(real*real + imag*imag)) / N);
>
> but how can I find the phase?
> (http://en.wikipedia.org/wiki/Frequency_spectrum say that it is possible
> but not how)

atan2(imag, real)

> If I have a silence period of 0.2 seconds before starting a sinus of 120
> Hz with an amplitude of 5000, will the silence be erased from the rebuilt
> file?
>
> Actually, in order to rebuild the original signal, I need to know the
> frequency, the amplitude, the start time and the end time of each sinus?

If you just perform the DFT and then the inverse DFT (properly
normalized; see above), then your output will be the same as the input
(modulo floating-point errors).

It sounds like you need to read some basic introductions to signal
processing and discrete Fourier transforms.  I'm sure you can find
several books on the subject in your library.  We also have some links
at www.fftw.org/links

Cordially,
Steven G. Johnson