Fast Fourier transform FFT? One line of python code is enough

Posted by Tremour on Sat, 23 Oct 2021 08:45:13 +0200

This article will introduce how to use one line of python code to realize the definition of base two time sampling FFT function.

Before we get to the point, let's warm up and use python to realize one line quick sorting. This is relatively easy. List derivation is a very convenient thing, so we just need:

quick_sort = lambda x :quick_sort([i for i in x if i<x[0]])+[i for i in x if i==x[0]]+quick_sort([i for i in x if i>x[0]]) if len(x) else []

It's very simple. You can do it with a recursion. Well, let's try to implement the FFT algorithm with one line of code.

Although it is one line, for the sake of the simplicity of the code, we first define the two mathematical limits of π and e (landlord's title party!), but after all, it is only a constant assignment, and only for the sake of the simplicity of the code, so please understand.  

import math
e = math.exp(1)
pi = math.pi

  Firstly, the first problem we encounter is that the basic two time FFT algorithm requires that the signal length must be the n-power of 2, otherwise we have to fill zero after the signal until the signal length becomes the minimum n-power greater than its original length. So how do you implement this function with only one line?

If we want to judge whether a positive integer is to the nth power of 2, first we may think of dividing it by 2 repeatedly until there is a remainder or it is equal to 1. But this is obviously not what one line of code can do. At this time, we can think of python's bit operation. As we all know, numbers are actually saved in binary in python, so we can use bit operation to judge whether a number is the N-power of 2. Here is a direct conclusion. If X & (x-1) is 0, it must be the n-th power of 2 (it should be noted that it is also satisfied when x=0, so it should be judged separately).

So next, how should we judge the n-th power greater than the minimum 2 of a number? The bin function can convert a number into a string of its binary number (format "0bxxxxx"). Next, you only need to judge the length of the string to know how much it is. Here is the written result directly.

y = lambda x:x + [0] * (2 ** (len(bin(len(x))) - 2) - len(x)) if (len(x) & (len(x) - 1)) and (len(x) > 2) else x

In this way, we get the function of automatically adding x to 0, but it is still thousands of miles away from the FFT algorithm. Therefore, we have to consider how to recursively implement FFT when the signal length meets the requirements.

FFT = lambda x :FFT(x + [0] * (2 ** (len(bin(len(x))) - 2) - len(x))) if (len(x) & (len(x) - 1)) and (len(x) > 2) else

First, judge that if the signal length is greater than 2 and the length does not meet the requirements, it will automatically fill 0 and call FFT recursively. Then, after this else, we can directly write our signal processing.

FFT = lambda x :FFT(x + [0] * (2 ** (len(bin(len(x))) - 2) - len(x))) if (len(x) & (len(x) - 1)) and (len(x) > 2) else [x1 + x2 * (len(x) / 2 - j - 0.5)/abs(len(x) / 2 - j - 0.5) for j ,[x1 ,x2] in enumerate([[x1 ,x2 * (e ** ((-1j) * 2 * pi * i / len(x)))] for i ,(x1 ,x2) in enumerate(zip(FFT([x[i] for i in range(len(x)) if not(i%2)]),FFT([x[i] for i in range(len(x)) if i%2])))] * 2)] if len(x)>2 else

The following long code may be a little difficult to understand, but it is also the core part. However, we analyze the part behind the first else separately. We know that the basic two-time sampling algorithm is an algorithm that converts the discrete Fourier transform of long signal into the discrete Fourier transform of short signal, and recurses until the signal length is 2.

Therefore, we judge that when the signal length is 2, we divide the even sequence and odd sequence of the signal, and then call it recursively again. In order to prevent the recursive operation of each layer from repeatedly calculating the next layer, resulting in the exponential growth of time complexity, we only calculate the results once and the number of times to be used after replication (that is, the role of * 2 there), It is packaged into a list by combining the enumerate and zip functions, so that we can use the list derivation we used before to deal with the Fourier transform results of the parity sequence returned by the next layer. Mention the role of the middle (len(x) / 2 - j - 0.5)/abs(len(x) / 2 - j - 0.5), which is responsible for separating the first half from the second half. When the i index is the first half of the sequence, its value is + 1, and vice versa. In this way, it replaces the if statement to realize the separate processing of the first and second half of the sequence. (if they are processed separately and connected with a + sign, the next layer FFT will have to be calculated twice, resulting in a very large increase in time complexity)

FFT = lambda x :FFT(x + [0] * (2 ** (len(bin(len(x))) - 2) - len(x))) if (len(x) & (len(x) - 1)) and (len(x) > 2) else [x1 + x2 * (len(x) / 2 - j - 0.5)/abs(len(x) / 2 - j - 0.5) for j ,[x1 ,x2] in enumerate([[x1 ,x2 * (e ** ((-1j) * 2 * pi * i / len(x)))] for i ,(x1 ,x2) in enumerate(zip(FFT([x[i] for i in range(len(x)) if not(i%2)]),FFT([x[i] for i in range(len(x)) if i%2])))] * 2)] if len(x)>2 else [x[0] + x[1],x[0] - x[1]] if len(x)==2 else x

The last step is to achieve success. When we recurse layer by layer until the signal length is 2, we can finally return its value directly! In addition, for the robustness of the code, when the signal length of x is 1, we also return its DFT result - itself. In this way, we have completed the FFT algorithm implemented with one line of code!

Use matalb to verify:

X = [1,4,4,1,4,1,1,4]
print(FFT(X))

-------------

The result is:
[(20+0j), (1.2426406871192857-3j), 0j, (-7.242640687119286+2.9999999999999996j), 0j, (-7.242640687119286-3j), 0j, (1.2426406871192857+3.0000000000000004j)]

  perfect!

Topics: Python Back-end