О классификации методов преобразования Фурье на примерах их программной реализации средствами Python26.09.2017 10:18

Введение


Публикации по методу Фурье условно можно разделить на две группы. Первая группа так называемых познавательных публикаций, например, [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()


Результат


27cc83b3f3b249338e91e7d708ce045b.JPG

4d27941bab774728a3b4908e61cd9d0d.JPG

Спектральный анализ периодических функций заключается в нахождении амплитуды А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


Результат


a490ac6dd2894def912796a4a7c8f97f.JPG

dcfc0d35cdb449a883cf9685581ed7be.JPG

ed033b61d7004e8a900d0450b385aa0e.JPG

Программная реализация спектрального анализа и синтеза с использованием функций 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


2f8b18f5f606405daacb0103143495a2.JPG

fb61c770a7c540fc892d1b4b286d5fc8.JPG

Фильтрация аналоговых сигналов


Под фильтрацией подразумевается выделение полезного сигнала из его смеси с мешающим сигналом — шумом. Наиболее распространенный тип фильтрации — частотная фильтрация. Если известна область частот, занимаемых полезным сигналом, достаточно выделить эту область и подавить те области, которые заняты шумом.

Программная реализация иллюстрирует технику фильтрации с применением БПФ. Сначала синтезируется исходный сигнал, представленный 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()


Результат


a8dcccb9b60e43f29d9c329f85741fc4.JPG

56a545bb8fd54e3b8fa43caed40e5942.JPG

Решение дифференциальных уравнений в частных производных


Алгоритм решения дифференциальных уравнений математической физики с использованием прямого и обратного БПФ приведен в [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()


Результат


3b5ac64fae48413b9db713b12b1bc674.JPG

Вывод


В статье приведена попытка классификации по областям применения методов преобразования Фурье.

Ссылки


  1. Математика на пальцах: давайте посчитаем хотя бы один ряд Фурье в уме.
  2. Простыми словами о преобразовании Фурье.
  3. Спектральный анализ сигналов нелинейных звеньев АСУ на Python.
  4. Математика в Python: Преобразование Фурье.
  5. Спектральный метод на примере простых задач матфизики.

© Habrahabr.ru