Transformada rapida de Fourier con FFTW

Autor: Mariano Maccarrone Aqui dejo un articulo que hace bastante queria escribir pero por una razon u otra no podia. Este articulo trata sobre la fft, fast fourier transform, o transformada rapida de fourier con la libreria fftw en C.

Instalacion

La instalacion es muy simple, en FreeBSD se puede instalar con sysinstal, si no compilarla es muy simple, se baja un archivo tar.gz se descomprime y con los clasicos comandos, ./configure, make y make install queda funcionando.

Uso

Este algoritmo es muy util para procesamiento de imagenes y filtrado digital pero en este articulo solo les voy a mostrar como usarla para hacer la transformada de un vector de una dimension. El ejemplo calcula un ejercicio sacado del libro "Matematicas avanzadas para ingenieria" de Glyn James que dice: Hallar usando fft la transformada de la sucesion g = {1,2,1,0}. Desde ya les adelanto que el resultado es {4,-2j,0,2j}. Seguido a esto calcularemos la transformada inversa a ver si nos da lo mismo.
Aca esta el programa:
#include <stdio.h>
#include <stdlib.h>
#include <fftw3.h>


main()
{
        int aux;
        double vec1[]={1,2,1,0};
        double vec2[]={0,0,0,0};
        double N=4;

        fftw_complex *in, *out, *in_inv, *out_inv;
        fftw_plan p;

        in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
        out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

        for (aux=0;aux<(int)N;aux++)
        {
                in[aux][0] = vec1[aux];
                in[aux][1] = vec2[aux];
        }

        p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

        fftw_execute(p);

        printf("Directa\n");
        for (aux=0;aux<(int)4;aux++)
        {
                printf("(%lf, %lf j)\n",out[aux][0],out[aux][1]);
        }

        in_inv = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
        out_inv = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

        for (aux=0;aux<(int)N;aux++)
        {
                in_inv[aux][0] = out[aux][0];
                in_inv[aux][1] = out[aux][1];
        }

        p = fftw_plan_dft_1d(N, in_inv, out_inv, FFTW_BACKWARD, FFTW_ESTIMATE);
        fftw_execute(p);

        printf("Inversa\n");
        for (aux=0;aux<(int)N;aux++)
        {
                printf("%lf,  %lf j)\n",out_inv[aux][0]/N,out_inv[aux][1]/N);
        }

        fftw_destroy_plan(p);
        fftw_free(in); fftw_free(out);
}
Comentarios del programa:

Los vectores vec1 y vec2 guardan la parte real y la imaginaria respectivamente de el vector a transformar. El tipo fftw_complex es simplemente esto
typedef double fftw_complex[2] 
donde el valor 1 corresponde a la parte real y el valor 2 a la imaginaria. La funcion fftw_plan_dft_1d crea la transformada, hay que pasarle el valor N, el puntero de entrada, el de salida, si es la fft directa se le pasa FFT_FORWARD y si es la inversa FFT_BACKWARD y FFT_ESTIMATE es un valor para optimizar el procesamiento de la transformada (ver manual en la pagina). La fft se ejecuta con la funcion fftw_execute y devuelve el resultado en el vector out.
Este programa lo compilamos con:
$ cc -lfftw3 -lm -o fft2 fft2.c
y si lo ejecutamos:
$ ./fft2
Directa
(4.000000, 0.000000 j)
(0.000000, -2.000000 j)
(0.000000, 0.000000 j)
(0.000000, 2.000000 j)
Inversa
1.000000,  0.000000 j)
2.000000,  0.000000 j)
1.000000,  0.000000 j)
0.000000,  0.000000 j)
v Como vemos el resultado es correcto y la mnversa tambien es correcta.
Por hoy me despido, para los que no sepan que es la fft tratare de escribir mas articulos explicando sobre su teoria y aplicaciones.

taxon:

Comentarios