О классификации методов преобразования Фурье на примерах их программной реализации средствами Python
Публикации по методу Фурье условно можно разделить на две группы. Первая группа так называемых познавательных публикаций, например, [1,2].
Вторая группа публикаций касается применения преобразований Фурье в технике, например, при спектральном анализе [3,4].
Ни в коем случае не умоляя достоинства этих групп публикации стоит признать, что без классификации, или хотя бы попытки осуществить такую классификацию, получить системное представление о методе Фурье, по моему мнению, затруднительно.
Задачи публикации
Провести классификацию методов преобразования Фурье на примерах их программной реализации средствами Python. При этом для облегчения чтения использовать формулы только в программном коде с соответствующими пояснениями.
Гармонический анализ и синтез
Гармоническим анализом называют разложение функции f(t), заданной на отрезке [0, Т] в ряд Фурье или в вычислении коэффициентов Фурье по формулам.
Гармоническим синтезом называют получение колебаний сложной формы путем суммирования их гармонических составляющих (гармоник).
#!/usr/bin/python # -*- coding: utf-8 -* from scipy.integrate import quad # модуль для интегрирования import matplotlib.pyplot as plt # модуль для графиков import numpy as np # модуль для операций со списками и массивами T=np.pi; w=2*np.pi/T# период и круговая частота def func(t):# анализируемая функция if t=pi f(t)=-cos(t)") plt.plot(q, F1, label='1 гармоника') plt.plot(q, F2 , label='2 гармоника') plt.plot(q, F3, label='3 гармоника') plt.xlabel("Время t") plt.ylabel("Амплитуда А") plt.legend(loc='best') plt.grid(True) F=np.array(a[0]/2)+np.array([0*t for t in q-1])# подготовка массива для анализа с a[0]/2 for k in np.arange(1,c,1): F=F+np.array([a[k]*np.cos(w*k*t)+b[k]*np.sin(w*k*t) for t in q])# вычисление членов ряда Фурье plt.figure() P=[func(t) for t in q] plt.title("Классический гармонический синтез") plt.plot(q, P, label='f(t)') plt.plot(q, F, label='F(t)') plt.xlabel("Время t") plt.ylabel("f(t),F(t)") plt.legend(loc='best') plt.grid(True) plt.show() Результат
Спектральный анализ периодических функций заключается в нахождении амплитуды Аk и фазы j k гармоник (косинусоид) ряда Фурье. Задача, обратная спектральному анализу, называется спектральным синтезом.
Программная реализация для спектрального анализа и синтеза без специальных функций NumPy для Фурье преобразования
#!/usr/bin/python # -*- coding: utf-8 -* from scipy.integrate import quad # модуль для интегрирования import matplotlib.pyplot as plt # модуль для графиков import numpy as np # модуль для операций со списками и массивами T=np.pi; w=2*np.pi/T# период и круговая частота def func(t):# анализируемая функция if t Результат
Программная реализация спектрального анализа и синтеза с использованием функций NumPy прямого быстрого преобразования Фурье – rfft и обратного преобразования – irfft
#!/usr/bin/python # -*- coding: utf-8 -* import matplotlib.pyplot as plt import numpy as np import numpy.fft T=np.pi;z=T/16; m=[k*z for k in np.arange(0,16,1)];arg=[];q=[]# 16 отсчётов на период в пи def f(t):# анализируемая функция if t 
Фильтрация аналоговых сигналов
Под фильтрацией подразумевается выделение полезного сигнала из его смеси с мешающим сигналом — шумом. Наиболее распространенный тип фильтрации — частотная фильтрация. Если известна область частот, занимаемых полезным сигналом, достаточно выделить эту область и подавить те области, которые заняты шумом.
Программная реализация иллюстрирует технику фильтрации с применением БПФ. Сначала синтезируется исходный сигнал, представленный 128 отсчетами вектора v. Затем к этому сигналу присоединяется шум с помощью генератора случайных чисел (функция np. random.uniform(0,0.5)) и формируется вектор из 128 отсчетов зашумленного сигнала.
#!/usr/bin/python # -*- coding: utf-8 -* import matplotlib.pyplot as plt import numpy as np from numpy.fft import rfft, irfft from numpy.random import uniform k=np.arange(0,128,1) T=np.pi;z=T/128; m=[t*z for t in k]#задание для дискретизации функции на 128 отсчётов def f(t):#анализируемая функция if t=0: q=1 else: q=0 return q v=[f(t) for t in m]#дискретизация исходной функции vs= [f(t)+np.random.uniform(0,0.5) for t in m]# добавление шума plt.figure() plt.title("Фильтрация аналоговых сигналов \n Окно исходной и зашумленной функций") plt.plot(k,v, label='Окно исходной функции шириной pi') plt.plot(k,vs,label='Окно зашумленной функции шириной pi') plt.xlabel("Отсчёты -k") plt.ylabel("Амплитуда А") plt.legend(loc='best') plt.grid(True) al=2# степень фильтрации высших гармоник fs=np. fft.rfft(v)# переход из временной области в частотную с помощью БПФ g=[fs[j]*FH(abs(fs[j])-2) for j in np.arange(0,65,1)]# фильтрация высших гармоник h=np.fft.irfft(g) # возврат во временную область plt.figure() plt.title("Фильтрация аналоговых сигналов \n Результат фильтрации") plt.plot(k,v,label='Окно исходной функции шириной pi') plt.plot(k,h, label='Окно результата фильтрации шириной pi') plt.xlabel("Отсчёты -k") plt.ylabel("Амплитуда А") plt.legend(loc='best') plt.grid(True) plt.show() Результат
Решение дифференциальных уравнений в частных производных
Алгоритм решения дифференциальных уравнений математической физики с использованием прямого и обратного БПФ приведен в [5]. Воспользуемся приведенными данными для программной реализации на Python решения дифференциального уравнения распространения тепла в стержне с применением преобразования Фурье.
#!/usr/bin/env python #coding=utf8 import numpy as np from numpy.fft import fft, ifft # функции быстрого прямого и обратного преобразования Фурье import pylab# модуль построения поверхности from mpl_toolkits.mplot3d import Axes3D# модуль построения поверхности n=50# стержень длиной 2 пи разбивается на 50 точек times=10# количества итераций во времени h=0.01# шаг по времени x=[r*2*np.pi/n for r in np.arange(0,n)]# дискретизация х w= np.fft.fftfreq(n,0.02)# сдвиг нулевой частотной составляющей к центру спектра k2=[ r**2 for r in w]# коэффициенты разложения u0 =[2 +np.sin(i) + np.sin(2*i) for i in x]# дискретизация начального распределения температуры вдоль стержня u = np.zeros([times,n])# матрица нулей размерностью 10*50 u[0,:] = u0 # нудевые начальные условия uf =np.fft. fft(u0) # коэффициенты Фурье начальной функции for i in np.arange(1,times,1): uf=uf*[(1-h*k) for k in k2] #следующий временной шаг в пространстве Фурье u[i,:]=np.fft.ifft(uf).real #до конца физического пространства X = np.zeros([times,n])# подготовка данных координаты х для поверхности for i in np.arange(0,times,1): X[i:]=x T = np.zeros([times,n])# подготовка данных координаты t для поверхности for i in np.arange(0,times,1): T[i:]=[h*i for r in np.arange(0,n)] fig = pylab.figure() axes = Axes3D(fig) axes.plot_surface(X, T, u) pylab.show() Результат
Вывод
В статье приведена попытка классификации по областям применения методов преобразования Фурье.
Ссылки
numpy.fft.rfft#
Compute the one-dimensional discrete Fourier Transform for real input.
This function computes the one-dimensional n-point discrete Fourier Transform (DFT) of a real-valued array by means of an efficient algorithm called the Fast Fourier Transform (FFT).
Parameters : a array_like
n int, optional
Number of points along transformation axis in the input to use. If n is smaller than the length of the input, the input is cropped. If it is larger, the input is padded with zeros. If n is not given, the length of the input along the axis specified by axis is used.
axis int, optional
Axis over which to compute the FFT. If not given, the last axis is used.
Normalization mode (see numpy.fft ). Default is “backward”. Indicates which direction of the forward/backward pair of transforms is scaled and with what normalization factor.
New in version 1.20.0: The “backward”, “forward” values were added.
The truncated or zero-padded input, transformed along the axis indicated by axis, or the last one if axis is not specified. If n is even, the length of the transformed axis is (n/2)+1 . If n is odd, the length is (n+1)/2 .
If axis is not a valid axis of a.
For definition of the DFT and conventions used.
The one-dimensional FFT of general (complex) input.
The n-dimensional FFT of real input.
When the DFT is computed for purely real input, the output is Hermitian-symmetric, i.e. the negative frequency terms are just the complex conjugates of the corresponding positive-frequency terms, and the negative-frequency terms are therefore redundant. This function does not compute the negative frequency terms, and the length of the transformed axis of the output is therefore n//2 + 1 .
When A = rfft(a) and fs is the sampling frequency, A[0] contains the zero-frequency term 0*fs, which is real due to Hermitian symmetry.
If n is even, A[-1] contains the term representing both positive and negative Nyquist frequency (+fs/2 and -fs/2), and must also be purely real. If n is odd, there is no term at fs/2; A[-1] contains the largest positive frequency (fs/2*(n-1)/n), and is complex in the general case.
If the input a contains an imaginary part, it is silently discarded.
>>> np.fft.fft([0, 1, 0, 0]) array([ 1.+0.j, 0.-1.j, -1.+0.j, 0.+1.j]) # may vary >>> np.fft.rfft([0, 1, 0, 0]) array([ 1.+0.j, 0.-1.j, -1.+0.j]) # may vary
Notice how the final element of the fft output is the complex conjugate of the second element, for real input. For rfft , this symmetry is exploited to compute only the non-negative frequency terms.
numpy.fft.rfft#
Compute the one-dimensional discrete Fourier Transform for real input.
This function computes the one-dimensional n-point discrete Fourier Transform (DFT) of a real-valued array by means of an efficient algorithm called the Fast Fourier Transform (FFT).
Parameters : a array_like
n int, optional
Number of points along transformation axis in the input to use. If n is smaller than the length of the input, the input is cropped. If it is larger, the input is padded with zeros. If n is not given, the length of the input along the axis specified by axis is used.
axis int, optional
Axis over which to compute the FFT. If not given, the last axis is used.
Normalization mode (see numpy.fft ). Default is “backward”. Indicates which direction of the forward/backward pair of transforms is scaled and with what normalization factor.
New in version 1.20.0: The “backward”, “forward” values were added.
The truncated or zero-padded input, transformed along the axis indicated by axis, or the last one if axis is not specified. If n is even, the length of the transformed axis is (n/2)+1 . If n is odd, the length is (n+1)/2 .
If axis is not a valid axis of a.
For definition of the DFT and conventions used.
The one-dimensional FFT of general (complex) input.
The n-dimensional FFT of real input.
When the DFT is computed for purely real input, the output is Hermitian-symmetric, i.e. the negative frequency terms are just the complex conjugates of the corresponding positive-frequency terms, and the negative-frequency terms are therefore redundant. This function does not compute the negative frequency terms, and the length of the transformed axis of the output is therefore n//2 + 1 .
When A = rfft(a) and fs is the sampling frequency, A[0] contains the zero-frequency term 0*fs, which is real due to Hermitian symmetry.
If n is even, A[-1] contains the term representing both positive and negative Nyquist frequency (+fs/2 and -fs/2), and must also be purely real. If n is odd, there is no term at fs/2; A[-1] contains the largest positive frequency (fs/2*(n-1)/n), and is complex in the general case.
If the input a contains an imaginary part, it is silently discarded.
>>> np.fft.fft([0, 1, 0, 0]) array([ 1.+0.j, 0.-1.j, -1.+0.j, 0.+1.j]) # may vary >>> np.fft.rfft([0, 1, 0, 0]) array([ 1.+0.j, 0.-1.j, -1.+0.j]) # may vary
Notice how the final element of the fft output is the complex conjugate of the second element, for real input. For rfft , this symmetry is exploited to compute only the non-negative frequency terms.