Learning Dip 5th Day--FFT algorithm and implementation of C language

• Why FFT is required

The first question is why it is necessary to create an FFT, simply speaking, for speed. We acknowledge that DFT is very useful, but we find that his speed is not very fast, 1D DFT primitive algorithm time complexity is O (n^2), this can be observed by the formula, for 2D DFT its time complexity is O (n^4), this speed is really difficult to accept, that is to say, When you calculate a 1024x768 image, you will wait for the approximate ... Probably... I haven't tried, anyway, for a long time. Do not believe yourself to try. Therefore, it is imperative to find a fast method to calculate the FFT. The following is the DFT formulaCalculates a 4-point DFT. The calculation amount is as follows:

• How to get an FFT

by observing the DFT formula, we find that the DFT computes each x[u] and traverses the complete x[n] so, our idea is to be able to get x[u by other x[u+m] (M is an integer, can be negative), that is, to make full use of the computational structure between them to build the present results. , this method is easily represented as an iterative form. This paper describes the FFT of base 2, and the number of discrete signal x[n] is 2 k-th square, that is, if the complete discrete signal has n components {X1,X2,X3....XN} n=1<<k.

• Fundamentals of Mathematics

The mathematical basis of the FFT is actually:                                     the rotation factor has the following properties:These properties allow us to use the previous calculation results to complete the subsequent calculations.

• 2-fft

observation: A discrete signal with only two values, assuming n=2, the use of properties 2-symmetry can be

• 4-fft

as in the above 2 points, the promotion to 4 points, will be so, where the action in the box for the above 2-FFT:The final algorithm:

• 8-fft

don't say much nonsense, as above:The calculation process is:The result is:The complete calculation process is shown in the following motion diagram. CSDN will not move, the net disk to see themselves:4-fft.gif

• 2^n-fft

Similarly , the extension to the 2^n, can get similar results, but we found that in order to make the output result is a sequence of results, the order of input has undergone a series of adjustments, and each step of the power parameter of the WN is also a change, we must get its change law to better write programs:Observe the change pattern of the WN:The node distance (set to Z) is a continuous addition to the z power of the WN, starting with the 0 power of the WN. The interval is then a Z-formula, again from 0 powers to Z-times. Repeat the above process:For example, the second level, the interval z=2, the node is a solid black point, node 0, 1, do not operate, node 2*w8^0, node 3*w8^2, node 4,5 do not do the operation. The inherent law is whether node I is multiplied by the WN: if (i%z== odd) node i*wn^ (step);the number of increments per step is determined by the current calculation level step=1<< (k (Total series)-I (current series))

• Reorder Input Parameters

The row I represents the array subscript, the blue number represents the actual data: its inherent law is that the subscript is the even-numbered will be mapped to their own subscript 1/2, the following table is an odd number is mapped to the second half of the array (SIZE_N/2) (Subscript-1) 1/2, the arranged data into the front and back two parts, Recursive process until there are only two elements. is stopped. The process is as follows:============================================================================================================================================================================================================================================= ======================= =========================================================================================================== ======================= =========================================================================================================== =======================

• Observation algorithm

at this point, we have basically done all the features of the FFT butterfly algorithm, and the next step is to implement it using code.

• Implementation code

main.c//fourer1d////Created by Tony on 14/11/16.//Copyright (c) 2014 Tony. All rights reserved.//#include <stdio.h> #include <stdlib.h> #include <time.h> #include <math.h >//mac M_PI has a macro definition in math.h, so here we select the row macro definition #ifndef m_pi#define m_pi 3.14159265358979323846#endif#define SIZE 1024*16# Define Value_max 1000//////////////////////////////////////////////////////////////////////defines a complex structure body////////////    struct complex_{double real; Double Imagin;}; typedef struct COMPLEX_ complex;//////////////////////////////////////////////////////////////////////defines a complex number calculation, including multiplication, addition, subtraction///////////////////////////////////////////////////////////////////void Add_complex (Complex * Src1,    Complex *src2,complex *dst) {dst->imagin=src1->imagin+src2->imagin; Dst->real=src1->real+src2->real;}    void Sub_complex (Complex * Src1,complex *src2,complex *dst) {dst->imagin=src1->imagin-src2->imagin; Dst-&gT;real=src1->real-src2->real;}    void Multy_complex (Complex * Src1,complex *src2,complex *dst) {double r1=0.0,r2=0.0;    Double i1=0.0,i2=0.0;    r1=src1->real;    r2=src2->real;    i1=src1->imagin;    i2=src2->imagin;    dst->imagin=r1*i2+r2*i1; Dst->real=r1*r2-i1*i2;} In the FFT there is a WN N-sub-item, which will be used continuously in the iteration, see the algorithm description/////// void Getwn (Double n,double Size_n,complex * DST) {double    X=2.0*m_pi*n/size_n;    Dst->imagin=-sin (x); Dst->real=cos (x);} Randomly generate an input that shows that the data part has been commented out//commented out the display part is the data display, Can observe results///////////////////////////////////////////////////////////////////void setinput (double * Data,int N) {//    printf ("SetInput signal:\n");    Srand ((int) time (0));        for (int i=0;i<size;i++) {data[i]=rand ()%value_max;    printf ("%lf\n", Data[i]); }}//////////////////////The DFT function is defined as a simple DFT definition, with time complexity O (n^2), and/or a two-layer loop in the following function, with step 1 for each layer, Size is n, so O (n*n),//commented off the display part of the data display, you can observe the results///////////////////////////////////////////////////////////////////void    DFT (double * Src,complex * Dst,int size) {clock_t start,end;        Start=clock ();        for (int m=0;m<size;m++) {double real=0.0;        Double imagin=0.0;            for (int n=0;n<size;n++) {double x=m_pi*2*m*n;            Real+=src[n]*cos (x/size);                imagin+=src[n]* (-sin (x/size));        } Dst[m].imagin=imagin;       Dst[m].real=real;        /* if (imagin>=0.0) printf ("%lf+%lfj\n", Real,imagin);    else printf ("%lf%lfj\n", real,imagin); */} end=clock (); printf ("DFT Use time:%lf for DataSize of:%d\n", (double) (End-start)/clocks_per_sec,size);} Define the IDFT function, which is a simple idft definition, time complexity O (n^2),// There are two layers of loops in the following function, each with a step of 1,size n for each layer, so O (n*n),/////////void IDFT (Complex *src,complex *dst,int size) {clock_t    Start,end;    Start=clock ();        for (int m=0;m<size;m++) {double real=0.0;        Double imagin=0.0;            for (int n=0;n<size;n++) {double x=m_pi*2*m*n/size;            Real+=src[n].real*cos (x)-src[n].imagin*sin (x);                    Imagin+=src[n].real*sin (x) +src[n].imagin*cos (x);        } real/=size;        Imagin/=size;            if (dst!=null) {dst[m].real=real;        Dst[m].imagin=imagin;        } if (imagin>=0.0) printf ("%lf+%lfj\n", Real,imagin);    else printf ("%lf%lfj\n", Real,imagin);    } end=clock ();    printf ("IDFT Use time:%lfs for DataSize of:%d\n", (double) (End-start)/clocks_per_sec,size); }//////////////////////////////////////////////////////////////////////defines the initialization data for the FFT, because the FFT data is remapped and the recursive structure is//////// int Fft_remap (double * Src,int size_n) {if (size_n==1) return 0;    Double * temp= (double *) malloc (sizeof (double) *size_n);        for (int i=0;i<size_n;i++) if (i%2==0) temp[i/2]=src[i];    else temp[(size_n+i)/2]=src[i];    for (int i=0;i<size_n;i++) src[i]=temp[i];    Free (temp);    Fft_remap (SRC, SIZE_N/2);    Fft_remap (SRC+SIZE_N/2, SIZE_N/2); return 1;} Define the FFT, see the algorithm description, the annotated display part of the data display, you can observe the results////// void FFT (double * Src,complex * dst,int size_n) {Fft_   Remap (SRC, size_n);    for (int i=0;i<size_n;i++)//printf ("%lf\n", Src[i]); clock_t start,end;    Start=clock ();    int k=size_n;    int z=0;    while (k/=2) {z++;    } k=z;    if (size_n!= (1<<k)) exit (0);    Complex * src_com= (complex*) malloc (sizeof (Complex) *size_n);    if (src_com==null) exit (0);      for (int i=0;i<size_n;i++) {  Src_com[i].real=src[i];    src_com[i].imagin=0;        } for (int i=0;i<k;i++) {z=0;                for (int j=0;j<size_n;j++) {if ((j/(1<<i))%2==1) {Complex wn;                Getwn (z, Size_n, &wn);                Multy_complex (&src_com[j], &wn,&src_com[j]);                z+=1<< (k-i-1);                Complex temp;                int neighbour=j-(1<< (i));                Temp.real=src_com[neighbour].real;                Temp.imagin=src_com[neighbour].imagin;                Add_complex (&temp, &src_com[j], &src_com[neighbour]);            Sub_complex (&temp, &src_com[j], &src_com[j]);        } else z=0; }}/* for (int i=0;i<size_n;i++) if (src_com[i].imagin>=0.0) {printf ("%lf+%lfj\n", src_        Com[i].real,src_com[i].imagin); } else printf ("%lf%lfj\n", Src_com[i].real,src_com[i].imagin); */for (int i=0;i<size_n;i++) {dst[i].imagin=src_com[i].imagin;dst[i].real=src_com[i].real;}    End=clock ();    printf ("FFT Use time:%lfs for DataSize of:%d\n", (double) (End-start)/clocks_per_sec,size_n); }////////////////////////////////////////////////////////////////////int Main (int argc, const char * argv[]) {double I    Nput[size];    Complex Dst[size];    SetInput (input,size);    printf ("\ n");    DFT (input, DST, SIZE);    printf ("\ n");    FFT (input, DST, SIZE);    IDFT (DST, NULL, SIZE); GetChar ();}

• Test results:

The time is S,fft advantage obviously

Learning Dip 5th Day--FFT algorithm and implementation of C language