Ctrl+K

Быстрое преобразование Фурье (FFT)

Быстрое преобразование Фурье (FFT)#

При любом гуглинге встроенной в питоне реализации дискретного Фурье, узнаете о функции np.fft.fft. Вот только она работает не по тем формулам, которые мы описали выше. Это легко понять, сравнив время работы.

import numpy as np
def DFT_slow(x):
    """Compute the discrete Fourier Transform of the 1D array x"""
    x = np.asarray(x, dtype=float)

    N = x.shape[0]
    n = np.arange(N)
    k = n.reshape((N, 1))

    M = np.exp(-2j * np.pi * k * n / N)

    return np.dot(M, x)

x = np.random.random(1024)
np.allclose(DFT_slow(x), np.fft.fft(x))

True

%timeit DFT_slow(x)
%timeit np.fft.fft(x)

104 ms ± 11.5 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
10.6 µs ± 495 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

Отличие в скорости в 10000 раз… С чем это связано?

Оказывается, Fast Fourier Transform (FFT) использует определенные симметрии комплексных экспонент, засчёт чего не приходится делать так много вычислительных операций. Если бы FFT был открыт на пару лет раньше, «возможно» гонку вооружений США и СССР удалось бы остановить в зародыше (https://www.youtube.com/watch?v=eQlSvfUuQNs).

Для нас главное запомнить отличие в скорости:

  • Обычный DFT работает за \(\mathcal{O}\left[N^2\right]\)

  • FFT работает за \(\mathcal{O}[N \log N]\)

При больших \(N\) разница колоссальная.

Теоретические основы FFT*#

Fast Fourier Transform (FFT) часто также называют алгоритмом бабочки или алгоритмом Кули-Тюки. Он состоит из элементарных шагов - операций «бабочка».

У него две реализации - с прореживанием по времени и по частоте.

В случае прореживания по частоте шаг выглядит так

image.png

Формула для вычисления «Бабочки» с прореживанием по времени: $\( Y_1 = X_1 + X_2 \\ Y_2 = (X_1 - X_2) \cdot {W^k_N} \)$

Обозначения:

\(X_1\), \(X_2\) – исходные точки;

\(Y_1\), \(Y_2\) – точки результата,

\(W^k_N = exp(\frac{-j\cdot 2\cdot π \cdot k}{N} )\) – комплексный коэффициент.

А сама схема алгоритма так (светлый кружок означает операцию «бабочка»)

image.png

Алгоритм с прореживанием по времени устроен похожим образом

image.png

Формула для вычисления «Бабочки» с прореживанием по времени: $\( Y_1 = X_1 + X_2 \cdot {W^k_N} \\ Y_2 = X_1 - X_2 \cdot {W^k_N} \)$

Обозначения:

\(X_1\), \(X_2\) – исходные точки;

\(Y_1\), \(Y_2\) – точки результата,

\(W^k_N = exp(\frac{-j\cdot 2\cdot π \cdot k}{N} )\) – комплексный коэффициент.

Простой пример FFT в python#

Сложим два синуса с разными частотами. Хотим восстановить частоты исходных синусов, зная только суммарную фуцнкцию.

# Python example - Fourier transform using numpy.fft method

import numpy as np
import matplotlib.pyplot as plt

# How many time points are needed i,e., Sampling Frequency
samplingFrequency   = 100;

# At what intervals time points are sampled
samplingInterval       = 1 / samplingFrequency;

# Begin time period of the signals
beginTime           = 0;

# End time period of the signals
endTime             = 10;

# Frequency of the signals
signal1Frequency     = 4;
signal2Frequency     = 7;

# Time points
time        = np.arange(beginTime, endTime, samplingInterval);

# Create two sine waves
amplitude1 = np.sin(2*np.pi*signal1Frequency*time)
amplitude2 = np.sin(2*np.pi*signal2Frequency*time)

# Create subplot
figure, axis = plt.subplots(4, 1)
figure.set_size_inches(8, 8)
plt.subplots_adjust(hspace=1)

# Time domain representation for sine wave 1
axis[0].set_title('Sine wave with a frequency of 4 Hz')
axis[0].plot(time, amplitude1)
axis[0].set_xlabel('Time')
axis[0].set_ylabel('Amplitude')

# Time domain representation for sine wave 2
axis[1].set_title('Sine wave with a frequency of 7 Hz')
axis[1].plot(time, amplitude2)
axis[1].set_xlabel('Time')
axis[1].set_ylabel('Amplitude')

# Add the sine waves
amplitude = amplitude1 + amplitude2

# Time domain representation of the resultant sine wave
axis[2].set_title('Sine wave with multiple frequencies')
axis[2].plot(time, amplitude)
axis[2].set_xlabel('Time')
axis[2].set_ylabel('Amplitude')

# Frequency domain representation
fourierTransform = np.fft.fft(amplitude)/len(amplitude)           # Normalize amplitude !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
fourierTransform = fourierTransform[range(int(len(amplitude)/2))] # Берём только половину спектра из-за симметрии!!!!!

#fourierTransform = fourierTransform[range(int(len(amplitude)/2))] # Exclude sampling frequency

tpCount     = len(amplitude)

values      = np.arange(int(tpCount/2))
#values      = np.arange(int(tpCount/2))

timePeriod  = tpCount/samplingFrequency

frequencies = values/timePeriod

# Frequency domain representation
axis[3].set_title('Fourier transform depicting the frequency components')
axis[3].plot(frequencies, abs(fourierTransform))
axis[3].set_xlabel('Frequency')
axis[3].set_ylabel('Amplitude')

plt.show()
../../_images/9dd91408cafd809e794df378dd46ba03f6b06fe817e7721998fc3f17d16b6649.png

Видим, что спектр симметричен относительно своей середины. Это всегда верно для действительного сигнала. Чтобы не считать зря «ненужную» вторую половину спектра, лучше использовать numpy.fft.rfft: https://numpy.org/doc/stable/reference/generated/numpy.fft.rfft.html

Содержание