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�e
#define SOUND_OFILE "nemo-gen.raw" // Fichier sortie
#define CHANNELS 1 // Canaux (2=st�r�o, 1=mono)
#define SAMPLE_BYTES 2 // Nombre d'octets par sample (2=16bits, 1=8bits)
#define N 44100 // Fr�quence d'�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�n�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�moire
printf("Allocating memory..\n");
in = fftw_malloc(sizeof(double) * N);
out = fftw_malloc(sizeof(fftw_complex) * N_OUT);
// Cr�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'�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�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�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�cup�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�quence = i , amplitude = ampli
ampli++;
printf("Freq= %u, ampli= %u\n", i, ampli);
}
i++;
}
// On lib�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