FFT (the most detailed and popular introductory manual)

FFT (the most detailed and popular introductory manual)


First of all, I need to declare that this article is slightly modified on the basis of reprinting. It was reprinted and modified with the permission of the original author ZLH_HHHH (Sori Hui elder sister) . Since I am a beginner and a scumbag in mathematics, so my elder sister Suggest that I write a handbook for the disabled, of course I readily accept it! ! !


The article is a bit urgent.
I hope to point out where there are mistakes, my learning FFT is a relatively slow process. Repeatedly during this period. I wrote this blog post only a very simple understanding. It can also help beginners when learning FFT. Somewhat biased. Avoid too much mental burden.

Go straight to the topic:

First of all: DFT itself is not responsible for multiplication between polynomials.

DFT is just a transformation.

FFT is a fast algorithm for DFT. (Divide and conquer improves efficiency)

<script id="MathJax-Element-4" type="math/tex; mode=display"></script>

Use FFT. We quickly transform the polynomial into a form that is convenient for calculation

Calculate the product of two polynomials in this convenient form of calculation.

At this time, although we have already got the target polynomial. But its form is not what we want

So then use the inverse operation of FFT to quickly transform back

<script id="MathJax-Element-7" type="math/tex; mode=display"></script>

Let us remember: The inverse operation of FFT is FFT 1

<script id="MathJax-Element-10" type="math/tex; mode=display"></script>

For a representation of n-degree polynomial A. The most common form (coefficient expression):

A(x)= i=0n 1aixi

The nth degree polynomial here refers to the polynomial with the highest degree term of n 1. 2

The complexity of polynomial multiplication in the form of elementary school students' hand-calculated coefficients is O(n2)

If we know n different points on the curve of a polynomial of degree n.

We can calculate this polynomial


Think of the coefficients as unknowns. List n equations. Solve these n coefficients.

In other words, given n points on the curve:

<(x0,y0),(x1,y1),(x2,y2),.....(xn 1,yn 1)><script id="MathJax-Element-18" type="math/tex; mode=display"><(x_0,y_0),(x_1,y_1),(x_2,y_2),.....(x_{n-1},y_{n-1})></script >

We can determine its coefficient form;

We call this method of expressing a polynomial: point value expression

He has many advantages. For example, a polynomial multiplication can be calculated in O(n) time

Give another polynomial:

B(x)= i=0n 1bixi

Let A(x) and B(x) take the same value of x to be expressed as:

<(x0,A0),(x1,A1),(x2,A2),...(x2n 2,A2n 2)><(x0,B0),(x1,B1),(x2, B2),.....(x2n 2,B2n 2)><script id="MathJax-Element-23" type="math/tex; mode=display"><(x_0,A_0),(x_1 ,A_1),(x_2,A_2),.....(x_{2n-2},A_{2n-2})>\<(x_0,B_0),(x_1,B_1),(x_2,B_2 ),.....(x_{2n-2},B_{2n-2})></script>

Here Ai=A(xi), Bi=B(xi)

Then A(x) B(x) is:

<(x0,A0B0),(x1,A1B1),(x2,A2B2),...(x2n 2,A2n 2B2n 2)><script id="MathJax-Element-27" type=" math/tex; mode=display"><(x_0,A_0B_0),(x_1,A_1B_1),(x_2,A_2B_2),.....(x_{2n-2},A_{2n-2}B_{2n -2})></script> 3

Very convenient. Go with the flow

The reason why A(x) and B(x) take more points. It is because the degree bound of A(x) B(x) increases.

Insufficient points will lead to less than 2n 1

For a polynomial of degree n. The complexity of the naive method of randomly finding n different points is O(n2)

<script id="MathJax-Element-34" type="math/tex; mode=display"></script>

Assume that n is an even number. Then we put A(x). Reorganized into two polynomials:

Let the subscript of its coefficient be a polynomial composed of even numbers as A[0](x):

A[0](x)=a0+a2x+a4x2+...+an 2xn2+1

Let the subscript of its coefficient be a polynomial composed of odd numbers as A[1](x):

A[1](x)=a1+a3x+a5x2+...+an 1xn2+1


A(xi)=A[0](x2i)+xiA[1](x2i) 4

It does not seem to optimize the calculation.

Because <x20,x21,....x2n 1><script id="MathJax-Element-41" type="math/tex; mode=display"> </script>

It is still possible to form n different numbers.

Then the scale we want to calculate will not be halved.

If we choose some x appropriately. Will the calculation be optimized?


<x0,x1,....xn2 1, x0, x1,.... xn2 1><script id="MathJax-Element-42" type="math/tex; mode=display" > </script>

After each number is squared. The resulting collection size is halved:

<x20,x21,....x2n2 1><script id="MathJax-Element-43" type="math/tex; mode=display"> </script>

Because ( x)2=x2

Then A( x)=A[0](x2) xA[1](x2)A(x)=A[0](x2)+xA[1](x2)

But this relationship will not pass on. The numbers obtained after squaring are all numbers not less than 0.

Recursion again returns to its original form. so awkward.

But this enlightens us. Take the opposite. Then lay it flat. The scale of the problem can be halved. But the recursion fails again. Is there a value that does not invalidate?

<script id="MathJax-Element-46" type="math/tex; mode=display"></script>

If for a set of numbers with an even number of elements. After each element is squared. The new collection obtained. After removing duplicate elements. The collection size can be halved. And if the new set is even. The new collection still satisfies the above properties. We call this set the nature of halving.

<script id="MathJax-Element-47" type="math/tex; mode=display"></script>

If we can quickly find a set of values of the independent variable x that satisfies the halving property.

Divide and conquer is feasible.

There is a kind of thing. Can meet our needs.

That is-n-th unit complex root:

(The use of complex numbers is really worrying. This is also my biggest mental burden when learning FFT)

We define i= 1

Then the shape is cos +i sin

The plural of has many good properties.

For example: (cos +i sin )k=cos k +i sin k

The above properties can be obtained by mathematical induction (relatively simple. No proof here

We remember:

Wn=cos2 n+i sin2 n

I will tell you when the polynomial takes the following value for x. Help us with DFT

<W0n,W1n,.....Wn 1n><script id="MathJax-Element-54" type="math/tex; mode=display"> </script>

The above n plurals are called. Complex root of n times unit. (The n-th unit complex root refers to these n numbers)

If we take Wkn as x. According to the decomposition just now:

A(Wkn)=A[0]( (Wkn)2)+WknA[1]((Wkn)2)

And (Wkn)2=W2kn=cos 2k 2 n+i sin2k 2 n=cos 2 kn2+i sin2 kn2=Wkn2

Because of the periodicity of trigonometric functions:

We have:

Wkn=cosk 2 n+i sink 2 n=cos(k mod n) 2 n+i sin(k mod n) 2 n=Wk mod nn

So: (Wkn)2=Wkn2=Wk mod n2n2

Then: <(W0n)2,(W1n)2,...(Wn 1n)2> is equivalent to <W0n2,W1n2,...Wn2 1n2,W0n2,W1n2,...Wn2 1n2>< script id="MathJax-Element-60" type="math/tex; mode=display"><(W_n^0)^2,(W_n^1)^2,...(W_n^{n-1} )^2>\is equivalent to/</script>

Surprised, not surprised, not surprised.

The value set of x takes the unit complex root. Not only satisfies the nature of halving. And there is a certain regularity. Keep consistent with the original collection.

This means that we only need to calculate:

<W0n2,W1n2,...Wn2 1n2>


The scope of the problem is halved. And if n is even. Squaring again still halves the size of the collection.

Note: A[0]k=A[0](Wkn2)=A[0]((Wkn)2)A[1]k=A[1](Wkn2)=A[1]((Wkn)2)

Then, when we get:

<A[0]0,A[0]1,..A[0]n2 1,A[1]0,A[1]1,...A[1]n2 1><script id= "MathJax-Element-65" type="math/tex; mode=display"> </script>

This means we got all A[0]((Win)2), A[1]((Win)2)


A( Wkn)=A[0]k modn2+WknA[1]k mod n2

You get: <A0,A1,A2,...An 1><script id="MathJax-Element-68" type="math/tex; mode=display"> </script>

Where Ak=A(Wkn)

of course. FFT calculation n should be an integer power of 2. This is because of the halving each time. We can fill the insufficient coefficients with 0.

The above process can be regarded as taking the unit complex root to calculate the target y vector, as shown in the following figure:

(W0n)0(W1n)0 (Wn 1n)0(W0n)1(W1n)1 (Wn 1n)1... ...(W0n)n 1(W1n) n 1 (Wn 1n)n 1 a0a1 an 1 1A1A

In the principle of as simple as possible. We are not talking about this matrix in particular.

Because I said that n equations about coefficients must be solvable. So the above matrix:

(W0n)0(W1n)0 (Wn 1n)0(W0n)1(W1n)1 (Wn 1n)1... ...(W0n)n 1(W1n) n 1 (Wn 1n)n 1

There must be an inverse matrix 5 . (Have a foundation in linear algebra, maybe no stranger)

Maybe you are still a little confused. It's ok:

What if we think of FFT as an algorithm for calculating the above special matrix multiplication.

Using FFT we can quickly get:

(W0n)0(W1n)0 (Wn 1n)0(W0n)1(W1n)1 (Wn 1n)1... ...(W0n)n 1(W1n) n 1 (Wn 1n)n 1 a0a1 an 1 (FFT) > 1A An 1

Simultaneously. I can tell you directly: D= (W0n)0(W1n)0 (Wn 1n)0(W0n)1(W1n)1 (Wn 1n)1... ...( W0n)n 1(W1n)n 1 (Wn 1n)n 1

The inverse matrix is:

E= 1n(W 0n)01n(W 1n)0 1n(W (n 1)n)01n(W 0n)11n(W 1n)1 1n(W (n 1)n)1... ...1n(W 0n)n 11n(W 1n)n 1 1n(W (n 1)n)n 1

Because I gave you the answer directly. . . So in doubt. We can calculate it.

Let us remember Y=E*D

Y[t][j]= k=0n 1E[t][k] D[k][j], where t,j [0,n 1]

Y[t][j]= k=0n 11nW tkn Wkjn=1n k=0n 1(Wj tn)k

The above can be regarded as the sum of geometric series.

When t==j: Y[t][j]=1n k=0n 11k=1

When t j, let u=j t:

Y[t][j]=1n k=0n 1(Wun)k=1n (Wun)n 1Wun 1 (geometric sequence formula)

Because: (Wun)n=Wunn=Wun mod nn=W0n=1

Therefore, when t j: Y[t][j]=0

and so. Y is the identity matrix. E and D are mutually inverse matrices.

and so:

(W0n)0(W 1n)0 (W (n 1)n)0(W0n)1(W 1n)1 (W (n 1)n) 1... ...(W0n)n 1(W 1n)n 1 (W (n 1)n)n 1 A0A1 An 1 (FFT ) > na0na1 nan 1

and so. FFT is actually obtained by changing Wkn of FFT to W kn.

And each element of the resulting target vector can be divided by n. 6.

Efficient implementation of FFT:

usually. We want the implementation of FFT to be as fast as possible. So FFT algorithms also use iterative structure instead of recursive structure.

The following describes a commonly used de-recursion method.

first of all. The FFT described above is a polynomial that can only handle the power of 2 integer powers.

So there is no special explanation. All have n=2k (If the coefficient is insufficient, use 0 to make up)

For the first step of recursion:

Input array <a0,a1,...an 1><script id="MathJax-Element-91" type="math/tex; mode=display"> </script>

Is decomposed into two parts, even and odd

<a0,a2,a4,...an 2>, <a1,a3,a5,...an 1><script id="MathJax-Element-92" type="math/tex; mode=display ">///,///</script>

The first character from the right in the subscript binary form is 0 and is assigned to the left set

The first character from the right in the subscript binary form is 1 and is assigned to the right set

More general. For the kth recursion:

The kth character from the right of the subscript is 0 and is assigned to the left set of the problem it is in

The kth character from the right of the subscript is 1 and is assigned to the right set of the problem it is in

Then for a number B in binary form. Express it as (note: n=2k ):

B= t=0k 1bt 2t

According to the previous analysis. The position that this number will be assigned to after recursion is:

P(B)= t=0k 1bt 2k 1 t, (because half of the position is deleted each time)

and so. The first k bits of the binary form of a number are symmetrical. It is the position after recursion.

For example: k=4, (0011)2 after symmetry >(1100)2

We perform the above retake of the array before calling the FFT.

Iteratively implement FFT from bottom to top


The FFT of the recursive structure is better understood. After understanding the FFT of the recursive structure.

It is not difficult to write a non-recursive structure

So here we make a summary of the FFT basic framework of the recursive structure.

Because the plural is used. So inevitable. Our following calculations are all in the range of complex numbers

For convenience. We use the symbol: Co(a,b) to represent the complex number a+bi

Then the polynomial:

A(x)= t=0n 1atxt= t=0n 1Co(at,0)xt

Wn=Co(cos2 n,sin2 n)

Wkn=Co(cos2k n,sin2k n)

Let: Y[t]=Co(at,0), t [0,n 1]

When n=1:


When n>1:

Let: Y[0][t]=Y[2t], Y[1][t]=Y[2t+1]

Calculate FFT(Y[0],n2) and FFT(Y[1],n2)

After the calculation is completed, let:

Y[t]=Y[0][t]+WtnY[1][t]Y[t+n2]=Y[0][t] WtnY[1][t] (where t [0,n2 1])

Return Y[];

Let me explain the above thought

The value of A at <W0n,W1n,...Wn 1n> <script id="MathJax-Element-111" type="math/tex"> </script> is stored in Y.

And return to the upper level. So for the current call:



When t [0,n2 1]:


A(Wt+n2n)=Y[t+n2]=A[0](Wt+n2n2)+Wt+n2nA[1](Wt+n2n2)=A[0](Wtn2) WtnA[1](Wtn2 )=Y[0][t] WtnY[1][t]


Give an iterative structure of the FFT template:

Define the plural:

struct Complex { double x,y; Complex ( double x1= 0.0 , double y1= 0.0 ) { x=x1; y=y1; } Complex operator -( const Complex &b) const { return Complex (xb.x,yb.y); } Complex operator +( const Complex &b) const { return Complex (x+bx,y+by); } Complex operator *( const Complex &b) const { return Complex (x*bx-y*by,x*b.y+y*bx); } }; Copy code

To recurse:

void change (Complex y[], int len) { int i,j,k; for (i= 1 ,j=len/ 2 ;i<len -1 ;i++) { if (i<j) swap (y[i],y[j]); k=len/2 ; while (j>=k) { j-=k; k/= 2 ; } j+=k; } } Copy code

Iterative structure of FFT

When on== 1: FFT
When on==-1: FFT
void FFT (Complex y[], int len, int on) { change (y,len); for ( int h = 2 ;h<=len;h<<= 1 ) { Complex wn (cos(on* 2 *pi/h),sin(on* 2 *pi/h)) ; for ( int j = 0 ;j<=len;j+=h) { Complex w ( 1 , 0 ) ; for ( int k=j;k<j+h/2 ;k++) { Complex u=y[k]; Complex t=w*y[k+h/2 ]; y[k]=u+t; y[k+h/2 ]=ut; w=w*wn; } } } if (on== -1 ) for ( int i= 0 ;i<len;i++) y[i].x/=len; } Copy code

Handicapped Handbook

"FFT" template
"Polynomial and Fast Fourier Transform" original text

  1. The naming part does not need to be too tangled, another more general term is to call the two processes forward and inverse of FFT as DFT and IDFT respectively.
  2. The principle of "Introduction to Algorithms" is used here, so we will redefine it.
  3. C(x)=A(x) B(x)=AiBi
  4. Because the order is reduced during the split, we need to substitute x2i and add xi before the odd number.
  5. The inverse matrix is the core of FFT . If there is no inverse matrix, the operation cannot be reversed.
  6. In the inverse operation, because both sides are multiplied by n for the convenience of the operation, the final result is n times the result we need, just divide by.