FFT parallel algorithm and its application-based on MPI (II.)

Two-dimensional serial FFT

Two-dimensional FFT is the implementation of the first line of the FFT to put the results back to the row, and then the column to do the FFT results in the column, after calculating all the columns, the result is the response of the two-dimensional FFT.

Code: main.c

#include "fft12_ifft12.h"
#include <stdio.h>
void Main ()
{   int i,j;
    COMPLEX td[16]={{1,4},{2,3},{3,0},{4,0},{5,0},{6,0},{7,0},{8,0},{9,0},{10,0},{11,0},{12,0},{13,0},{14,0},{15,0} , {16,0}};
    COMPLEX fd[16];
    FFT2 (TD, 4, 4,FD);
    for (i=0;i<4;i++)
    {for   (j=0;j<4;j++)
        {   printf ("%10.6f+j%10.6f,  ", FD[I*4+J].RE,FD [i*4+j].im);
        }
        printf ("\ n");}
}

Complex_oper.h

#ifndef com_h_included
#define com_h_included
#include <math.h>
#include <malloc.h>
# Include <string.h>

#define PI (double) 3.14159265359/

* plural definition */
typedef struct
{
    double re;
    double im;
} COMPLEX;


/* Complex plus operation *
/static COMPLEX Add (COMPLEX C1, COMPLEX C2)
{
    COMPLEX C;
    C.re=c1.re+c2.re;
    c.im=c1.im+c2.im;
    return c;
}


/* Complex minus operation *

/static COMPLEX Sub (COMPLEX C1, COMPLEX C2)
{
    COMPLEX C;
    C.re=c1.re-c2.re;
    c.im=c1.im-c2.im;
    return c;
}

/* Complex multiplication operation *

/static COMPLEX Mul (COMPLEX C1, COMPLEX C2)
{
    COMPLEX C;
    c.re=c1.re*c2.re-c1.im*c2.im;
    c.im=c1.re*c2.im+c2.re*c1.im;
    return c;
}

#endif

Fft12_ifft12.cpp

#include <math.h> #include <malloc.h> #include <string.h> #include "fft12_ifft12.h"/* fast-pay-in-a-TD-last-field value,
    FD is the frequency domain value, power is 2 power number */void FFT (COMPLEX * TD, COMPLEX * FD, int power) {int count;
    int i,j,k,bfsize,p;
    Double angle;

    COMPLEX *w,*x1,*x2,*x;
    /* Calculate the pay-in-sorrow transform points */count=1<<power;
    /* Allocate required memory for operation */w= (COMPLEX *) malloc (sizeof (COMPLEX) *COUNT/2);
    x1= (COMPLEX *) malloc (sizeof (COMPLEX) *count);
    X2= (COMPLEX *) malloc (sizeof (COMPLEX) *count);
        /* Calculate weighted coefficients */for (i=0;i<count/2;i++) {angle=-i*pi*2/count;
        W[i].re=cos (angle);
    W[i].im=sin (angle);
    }/* Writes the time domain point to the memory */memcpy (x1,td,sizeof (COMPLEX) *count); /* Butterfly operation */for (k=0;k<power;k++) {j=0;j<1<<k;j++) {bfsize=1<< (POW
            ER-K);
                for (i=0;i<bfsize/2;i++) {p=j*bfsize;
                X2[i+p]=add (X1[I+P],X1[I+P+BFSIZE/2]); X2[i+p+bfsize/2]=mul (Sub (X1[I+P],X1[I+P+BFSIZE/2]), w[i* (1<<k)]);
        }} x=x1;
        x1=x2;
    X2=x;
        }/* Reorder */for (j=0;j<count;j++) {p=0;
        for (i=0;i<power;i++) {if (j& (1<<i)) p+=1<< (power-i-1);
    } Fd[j]=x1[p];
    }/* Release memory */free (W);
    Free (X1);
Free (X2);
    }///Fast pay-in-sorrow inverse transform, using fast pay-in-mourning to transform FD into a frequency domain value, TD last field value, Power 2 power number */void IFFT (COMPLEX *fd, COMPLEX *td, int power) {int i,count;

    COMPLEX *x;
    /* Calculate the pay-in-sorrow inverse transform points */count=1<<power;
    /* Allocate required memory for operation */x= (COMPLEX *) malloc (sizeof (COMPLEX) *count);
    /* Write the frequency domain point to the memory */memcpy (x,fd,sizeof (COMPLEX) *count);
    /* Conjugate of frequency domain points */for (i=0;i<count;i++) {x[i].im=-x[i].im;
    }/* Invoke fast pay-in-mourning transform */FFT (x,td,power);
        /* Conjugate of the time domain point */for (i=0;i<count;i++) {td[i].re/=count;
    Td[i].im=-td[i].im/count;
}/* Release memory */free (x); }



/********************************************************* * Function Name: * Fourier () * * parameter: * complex* TD-input time domain sequence * Long lwidth-image Width *
 Long Lheight-Image Height * complex* FD-Output frequency domain sequence * return value: * BOOL-successfully Returns TRUE, otherwise returns false.
 * * Description: * This function is a two-dimensional fast Fourier transform. * ************************************************************************/void FFT2 (COMPLEX * TD, Long lwidth, long l
    Height, COMPLEX * FD) {COMPLEX *tempt, *tempf;
    Cyclic variable long i;
    Long J;


    Long K;
    The width and height of the Fourier transform (2 of the whole number of square) long w = 1;
    Long h = 1;
    int wp = 0;

    int hp = 0;
        Calculates the width and height of the Fourier transform (integer 2) while (W < lwidth) {W *= 2;
    wp++;
        } while (H < lheight) {h *= 2;
    hp++;
    }//Allocate memory tempt = (COMPLEX *) malloc (sizeof (COMPLEX) *h);

    Tempf = (COMPLEX *) malloc (sizeof (COMPLEX) *h); Fast Fourier transform for y direction//rgb/*for (i = 0; i < w * 3; i++) {//Extract data for (j = 0; J < H

        J + +) Tempt[j] = td[j * W * 3 + I];//RGB//one-dimensional fast Fourier transform FFT (tempt, Tempf, HP);
    Save transform results for (j = 0; J < H; j + +) Td[j * W * 3 + i] = tempf[j]; } *///grayscale for (i = 0; i < W; i++) {//Extract data for (j = 0; J < H; j + +) {Tem PT[J] = td[j * w + i];}

        RGB//One-dimensional fast Fourier transform FFT (tempt, Tempf, HP);
    Save transform results for (j = 0; J < H; j + +) {Td[j * w + i] = tempf[j];}
    }//Release memory free (tempt);

    Free (TEMPF);
    Allocate memory tempt = (COMPLEX *) malloc (sizeof (COMPLEX) *w);
    Tempf = (COMPLEX *) malloc (sizeof (COMPLEX) *w);
            Fast Fourier transform in x direction//RGB/* for (i = 0; i < H; i++) {for (k = 0; k < 3; k++) {

            Extract data for (j = 0; J < W; j + +) Tempt[j] = td[i * W * 3 + J * 3 + K];

            One-dimensional FFT with fast Fourier transform (tempt, Tempf, WP);
         Save Transform results   for (j = 0; J < W; j + +) Fd[i * W * 3 + J * 3 + K] = Tempf[j];
            }} *//Grayscale for (i = 0; i < H; i++) {//Extract data for (j = 0; J < W; J + +)

            {Tempt[j] = td[i * W + j];}

            One-dimensional FFT with fast Fourier transform (tempt, Tempf, WP);

    Save transform results for (j = 0; J < W; J + +) {Fd[i * w + j] = Tempf[j];}
    }//Release memory free (tempt);


Free (TEMPF);   }/************************************************************************* * * Function Name: * Ifourier () * * Parameter: *  Lpbyte TD-Spatial data returned * Long lwidth-width of airspace image * Long lheight-height of airspace image * complex* FD-
 Input frequency domain Data * * return value: * BOOL-successfully Returns TRUE, otherwise returns false.
 * * Description: * This function is a two-dimensional fast Fourier inverse transformation. * ************************************************************************/void IFFT2 (COMPLEX *td, long lWidth, long
    Lheight, COMPLEX * FD) {COMPLEX *tempt, *tempf; //cyclic variable long i;
    Long J;



    Long K;
    The width and height of the Fourier transform (2 of the whole number of square) long w = 1;
    Long h = 1;
    int wp = 0;

    int hp = 0;
        Calculates the width and height of the Fourier transform (integer 2) while (W < lwidth) {W *= 2;
    wp++;
        } while (H < lheight) {h *= 2;
    hp++;

    }//Calculates the number of bytes per line of the image//long Llinebytes = widthbytes (Lwidth * 24);
    Allocate memory tempt = (COMPLEX *) malloc (sizeof (COMPLEX) *w);

    Tempf = (COMPLEX *) malloc (sizeof (COMPLEX) *w);
            Fast Fourier transform in x direction//RGB/* for (i = 0; i < H; i++) {for (k = 0; k < 3; k++) {

            Extract data for (j = 0; J < W; j + +) Tempf[j] = fd[i * W * 3 + J * 3 + K];

            One-dimensional fast Fourier transform IFFT (Tempf, Tempt, WP);
        Save transform results for (j = 0; J < W; j + +) Fd[i * W * 3 + J * 3 + K] = Tempt[j];
}} *//Grayscale for (i = 0; i < H; i++) {//Extract data            for (j = 0; J < W; j + +) Tempf[j] = fd[i * W + j];

            One-dimensional fast Fourier transform IFFT (Tempf, Tempt, WP);

    Save transform results for (j = 0; J < W; j + +) Fd[i * W + j] = Tempt[j];
    }//Release memory free (tempt);

    Free (TEMPF);
    Tempt = (COMPLEX *) malloc (sizeof (COMPLEX) *h);

    Tempf = (COMPLEX *) malloc (sizeof (COMPLEX) *h); Fast Fourier transform for y direction//rgb/* for (i = 0; i < w * 3; i++) {//Extract data for (j = 0; J < H; j

        + +) Tempf[j] = fd[j * W * 3 + i];

        One-dimensional fast Fourier transform IFFT (Tempf, tempt, HP);
    Save transform results for (j = 0; J < H; j + +) Td[j * W * 3 + i] = tempt[j]; } *///grayscale for (i = 0; i < W; i++) {//Extract data for (j = 0; J < H; j + +) Tem

        PF[J] = fd[j * w + i];

        One-dimensional fast Fourier transform IFFT (Tempf, tempt, HP); Save transform results for (j = 0; J < H; j + +) td[J * W + i] = tempt[j];
    }//Release memory free (tempt);

    Free (TEMPF); /*for (i = 0; i < H; i++) {for (j = 0; J < W * 3, J + +) {if (I < lheight &AMP;&A mp
        J < llinebytes) * (TD + i * llinebytes + j) = Fd[i * W * 3 + j].re;
 }
    }
*/

}

fft12_ifft12.h

#ifndef com_h_included
#define com_h_included

#define PI (double) 3.14159265359/

* plural definition */
typedef struct
{
    double re;
    double im;
} COMPLEX;


/* Complex plus operation *
/static COMPLEX Add (COMPLEX C1, COMPLEX C2)
{
    COMPLEX C;
    C.re=c1.re+c2.re;
    c.im=c1.im+c2.im;
    return c;
}


/* Complex minus operation *

/static COMPLEX Sub (COMPLEX C1, COMPLEX C2)
{
    COMPLEX C;
    C.re=c1.re-c2.re;
    c.im=c1.im-c2.im;
    return c;
}

/* Complex multiplication operation *

/static COMPLEX Mul (COMPLEX C1, COMPLEX C2)
{
    COMPLEX C;
    c.re=c1.re*c2.re-c1.im*c2.im;
    c.im=c1.re*c2.im+c2.re*c1.im;
    return c;
}
void FFT (COMPLEX * TD, COMPLEX * FD, int power);
void IFFT (COMPLEX *fd, COMPLEX *td, int power);
void FFT2 (COMPLEX * TD, Long lwidth, long lheight, COMPLEX  * FD);
void IFFT2 (COMPLEX *td, Long lwidth, long lheight, COMPLEX * FD);

#endif