The statement

First of all, I need to declare that this article is reproduced on the basis of a little modified, after the original author ZLH_HHHH (Zuolihui senior student) permission just reprinted and modified, because I am a beginner, and is a number of poor students I, so my senior student suggested that I write a disabled manual, I of course is readily accepted!!

Body:

The article was written in a little hurry. Any mistakes please point out that FFT is a slow process for me to learn. It goes back and forth. I am writing this blog post with a very superficial understanding. At the same time can also help beginners in learning FFT. A little too much. Avoid too much mental baggage.

Let’s get straight to the point:

First: the DFT itself is not responsible for multiplying polynomials.

The DFT is just a transformation.

FFT is the fast algorithm of DFT. (Divide and conquer improves efficiency)

By using FFT. We quickly transform the polynomial into a form that is computationally convenient

The product of two polynomials is computed in this convenient form.

At this point we have our target polynomial. But the form is not what we want

So then the inverse operation of FFT is used to quickly transform back

We remember that the inverse of FFT is FFT−1

For A representation of A to the NTH degree. The most common form (coefficient expression) :

A (x) = ∑ aixi I = 0 n – 1

NTH degree polynomials here refer to polynomials of n−1 highest degree2

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

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

We can figure out this polynomial

According to?

Think of coefficients as unknowns. Let’s write out n equations. Solve for these n coefficients.

That is, 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 – valued representation

He has many good qualities. For example, you can compute a polynomial multiplication in order n time

Give me another polynomial:

B (x) = ∑ I = 0 n – 1 bixi

Let the point value of A(x) and B(x) with the same x be expressed as:

<(x0,A0),(x1,A1),(x2,A2),….. , A2n (x2n – 2-2) > < (x0, B0), (x1, B1) and (x2, B2),… , B2n (x2n – 2-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>

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

A (x) ∗ B (x) for:

<(x0,A0B0),(x1,A1B1),(x2,A2B2),….. B2n, A2n (x2n – 2-2-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

It’s very convenient. Let nature take its course;

That’s why I’m taking more points for A of x and B of x. That’s because the bounds on A(x)∗B(x) increase.

If you don’t take enough points, you don’t get to 2n−1

For a polynomial of NTH degree. The naive method for randomly finding n different points is order n squared.

Let’s say n is even. So let’s take A of x. Recombination into two polynomials:

Let the polynomial A[0](x) composed of its coefficients with even subscripts be:

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

Let the polynomial A[1](x) composed of its coefficients with odd subscripts be:

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

So there are:

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

It doesn’t seem to optimize.

Because the < x20, x21,… X2n – 1 > < script id = “MathJax – Element – 41” type = “math/Tex; mode=display”> </script>

It’s still possible to form n different numbers.

So we’re not going to cut the size of the calculation in half.

If we pick some x’s properly. Is it going to optimize the operation.

Such as:

<x0,x1,…. Xn2-1 – x0, – x1,… – xn2-1 > < script id = “MathJax – Element – 42” type = “math/Tex; mode=display”> </script>

You square the numbers. The resulting set 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 the relationship doesn’t pass on. If you square it, you get all the numbers that are at least 0.

Again, the recursion is back to the original form. Very embarrassed.

But it enlightens us. Take the negative. Then lay it flat. You can cut the size of the problem in half. But once again the recursion doesn’t work. Is there a value that doesn’t expire.

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

If I have an even number of numbers. Each element squared. I get a new set. After removing duplicate elements. The collection size can be halved. And I get a new set if it’s even. The new set still satisfies these properties. Let’s say this set has the property of being cut in half.

If we can quickly find a set of values of x that satisfy the halved property.

Divide-and-conquer works.

There is a thing. It can meet our needs.

That’s — n times the complex root of the unit:

(Plural is a real concern. This is also my biggest mental burden when LEARNING FFT.

We define I =−1‾ √

So cos theta + I ∗sin theta

The complex number of is has many nice properties.

For example :(cosθ+ I ∗sinθ)k=cos kθ+ I ∗sin kθ

The above properties can be obtained by mathematical induction. I’m not going to prove it here

We remember:

Wn = cos2 PI n + I ∗ sin2 PI n

I’m going to tell you that the polynomial is at x values below. It helps us do the DFT

<W0n,W1n,….. Wn – 1 n > < script id = “MathJax – Element – 54” type = “math/Tex; mode=display”> </script>

The n complex numbers above are called. N times the unit complex root. (NTH unit complex root refers to this n number.)

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

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

And (Wkn) 2 = W2kn = cos 2 k ∗ 2 PI n + I ∗ sin2k ∗ 2 PI n = cos 2 PI kn2 + I ∗ sin2 PI kn2 = Wkn2

Because 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

Is: “(W0n) 2, (W1n) 2,… (Wn−1n)2> is equivalent to <W0n2,W1n2… Wn2-1 n2, W0n2, W1n2,… Wn2-1 n2 > < script id = “MathJax – Element – 60” type = “math/Tex; mode=display”><(W_n^0)^2,(W_n^1)^2,… (W_n^{n-1})^2> = \\ </script>

It’s not a surprise, it’s not a surprise.

The set of values of x takes the unit complex root. Not only does it satisfy the property of half. And with some regularity. Same as the original set.

This means that we just need to compute:

<W0n2,W1n2,… Wn2-1 n2 >

A [0], A [1]

The scope of the problem halved. And if n is even. Square again and the set size is halved.

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

So, when we get:

<A[0]0,A[0]1,.. A [0] n2-1, A [1] 0, A [1] 1,… A [1] (n – 1 > < script id = “MathJax – Element – 65” type = “math/Tex; mode=display”> </script>

That means we get all of A[0]((Win)2),A[1]((Win)2)

Is:

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

<A0,A1,A2… The An – 1 > < script id = “MathJax – Element – 68” type = “math/Tex; mode=display”> </script>

The Ak = A (Wkn)

Of course. FFT is evaluated as n to the power of 2. That’s because it’s halved each time. We can fill in the missing coefficients with zeros.

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

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ (W0n) 0 (W1n) 0 ⋮ (Wn – 1, n) 0 (W0n) 1 (W1n) 1 ⋮ (Wn – 1, n) 1…. ⋱. (W0n) n – 1 (W1n) n – 1 ⋮ (Wn – 1, n) n – 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ a0a1 ⋮ the an – 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ a0a1 ⋮ The An – 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥

Keep it as simple as possible. We’re not talking about this matrix in particular.

Because we said that n equations of coefficients must be solvable. So the matrix above:

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ (W0n) 0 (W1n) 0 ⋮ (Wn – 1, n) 0 (W0n) 1 (W1n) 1 ⋮ (Wn – 1, n) 1…. ⋱. (W0n) n – 1 (W1n) n – 1 ⋮ (Wn – 1, n) n – 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥

There must be an inverse matrix5. (If you have a background in linear algebra, maybe you’re familiar with it.)

Maybe you’re still a little confused. It doesn’t matter:

So if we think of FFT as an algorithm to compute that particular matrix multiplication above.

Using FFT, we can quickly get:

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ (W0n) 0 (W1n) 0 ⋮ (Wn – 1, n) 0 (W0n) 1 (W1n) 1 ⋮ (Wn – 1, n) 1…. ⋱. (W0n) n – 1 (W1n) n – 1 ⋮ (Wn – 1, n) n – 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ a0a1 ⋮ the an – 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ – (FFT) – > ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ A0A1 ⋮ the An – 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥

At the same time. I can tell you: D = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ (W0n) 0 (W1n) 0 ⋮ (Wn – 1, n) 0 (W0n) 1 (W1n) 1 ⋮ (Wn – 1, n) 1…. ⋱. (W0n) n – 1 (W1n) n – 1 ⋮ (Wn – 1, n) n – 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥

The inverse matrix of is:

E = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 1 n (W – 0 n) n (W – 1 n) 0 01 ⋮ 1 n (W – (n – 1) n) n (n) W – 0 01 11 n – 1 n (W) 1 ⋮ 1 n (W – (n – 1) n) 1…. ⋱. 1 n (W – 0 n) n – 11 n (W – 1, n) n – 1 ⋮ 1 n (W – (n – 1) n) n – 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

Because I gave you the straight answer… So if you doubt it. So we can calculate it.

We 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 = 0 n nw – 11 – TKN ∗ Wkjn = 1 n ∑ k = 0 n – 1 (Wj – tn) k

You can view this as a sum of geometric sequences.

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 (formula of geometric sequence)

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

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

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

So:

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ (W0n) 0 (W – 1 n) 0 ⋮ (W – n) (n – 1) 0 (W0n) 1 (W – 1 n) 1 ⋮ (W – (n – 1) n) 1…. ⋱. (W0n) n – 1 (W – 1 n) n – 1 ⋮ (W – (n – 1) n) n – 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ A0A1 ⋮ the An – 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ – (FFT) – – > ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ na0na1 ⋮ nan – 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥

So. FFT minus is just Wkn of FFT becomes W minus kn.

And you get the target vector and you just subtract n from each element6

Efficient implementation of FFT:

Usually. We want FFT to be implemented as quickly as possible. So FFT algorithms also use iterative structures rather than recursive structures.

The following is a common method of de-recursion.

In the first place. The FFT introduced above is a polynomial that can only handle the degree bounds of an integer power of 2

So not specifically. They all have n=2k.

For the first step of recursion:

Input array < A0, A1… The an – 1 > < script id = “MathJax – Element – 91” type = “math/Tex; mode=display”> </script>

It’s broken down into even and odd numbers

<a0,a2,a4,… The an – 2 >, < a1, a3, a5,… The an – 1 > < script id = “MathJax – Element – 92” type = “math/Tex; mode=display”> \ \ \ ,\ \ \ </script>

Subscript binary form the first character of 0 from the right is assigned to the left set

Subscript binary form the first character of 1 from the right is assigned to the right set

More general. For the KTH recurrence:

The KTH character zero from the subscript right is assigned to the left set of its problem

The KTH character from the subscript right, 1, is assigned to the right set of its problem

So for a base 2 number B. It is expressed as (note: n=2k) :

B = ∑ t = 0 k – 1 bt ∗ 2 t

According to the previous analysis. The number will be recursively assigned to:

P(B)=∑t= 0K −1bt∗ 2K −1−t, (because half the positions are deleted each time)

So. For the binary form of a number, the first k bits are symmetric. That’s the recursion position.

For example, k=4, (0011)2− Symmetric −>(1100)2

We do one of these rebeats on the array before calling FFT.

FFT can be implemented iteratively from the bottom up

Conclusion:

Because the recursive structure of FFT is easier to understand. And then we understand the recursive structure of FFT.

It’s not hard to write non-recursive structures

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

Because we’re using complex numbers. So it’s inevitable. We’re going to do everything in the complex number range

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

Then the polynomial:

A (x) = = 0 n – 1 atxt = ∑ ∑ t t = 0 n – 1 co xt (ats, 0)

Wn = Co (cos2 PI n, sin2 PI n)

Wkn = Co (cos2k PI n, sin2k PI n)

Make: [t] Y = Co (ats, 0), t ∈ (0, n – 1)

When n = 1:

FFT(Y[],n)=Y[0]

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 complete, make:

[t] Y = Y [0] [t] + WtnY [1] [t] [t + n2] Y = Y [0] [t] – WtnY [1] [t] t ∈ ([0, n2-1))

Returns the Y [];

— — — — — — — — — — — —

Let me explain this idea

A. in the < W0n, W1n,… Wn-1n > <script id=” mathJax-element-111 “type=”math/ Tex “> </script> the values at the point are stored in Y.

And returns to the previous level. So for the current call:

A[0](Wtn2)=Y[0][t]A[1](Wtn2)=Y[1][t]

A(Wtn)=A[0]((Wtn)2)+WtnA[1]((Wtn)2)=A[0](Wtn2)+WtnA[1](Wtn2)

When t ∈ [0, n2-1) :

A(Wtn)=Y[t]=A[0](Wtn2)+WtnA[1](Wtn2)=Y[0][t]+WtnY[1][t]

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

template

An FFT template for an iterative structure is given:

Defining complex numbers:

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(x-b.x,y-b.y);
    }
    Complex operator+ (const Complex &b)const
    {
        return Complex(x+b.x,y+b.y);
    }
    Complex operator* (const Complex &b)const
    {
        return Complex(x*b.x-y*b.y,x*b.y+y*b.x); }};Copy the code

Recursive:

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 the code

FFT of iterative structures

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(intk=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]=u-t; w=w*wn; }}}if(on==- 1)
        for(int i=0; i<len; i++) y[i].x/=len; }Copy the code

Manual for persons with Disabilities

“FFT” template “Polynomial and fast Fourier Transform” original text

There is no need to be too tangled in the naming part. Another common appellation is to call the two processes of FFT DFT and IDFT respectively. ↩
I’m going to use the same principle that I used in Introduction to Algorithms, so let me redefine it. ↩
C (x) = A (x) ∗ B (x) = AiBi ↩
Because we dropped it when we split it, we’re going to plug in x2i, and we’re going to add x sub I to the odd case. ↩
The inverse matrix is the core of FFT−, and the inverse operation cannot be performed without the inverse matrix. ↩
And in the reverse case, because we multiplied both sides by n just to make it easier, we’re going to get n times what we want, divided by that. ↩

FFT (the most detailed and accessible how-to manual)