О классификации методов преобразования Фурье на примерах их программной реализации средствами Python
https://habr.com/ru/articles/338704/https://habr.com/ru/articles/338704/
26KPython*Математика*Разработка под Windows*Туториал
Введение
Публикации по методу Фурье условно можно разделить на две группы. Первая группа так называемых познавательных публикаций, например, [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<np.pi:
p=np.cos(t)
else:
p=-np.cos(t)
return p
def func_1(t,k,w):# функция для расчёта коэффициента a[k]
if t<np.pi:
z=np.cos(t)*np.cos(w*k*t)
else:
z=-np.cos(t)*np.cos(w*k*t)
return z
def func_2(t,k,w):#функция для расчёта коэффициента b[k]
if t<np.pi:
y=np.cos(t)*np.sin(w*k*t)
else:
y=-np.cos(t)*np.sin(w*k*t)
return y
a=[];b=[];c=4;g=[];m=np.arange(0,c,1);q=np.arange(0,2*np.pi,0.01)# подготовка списков для численного анализа
a=[round(2*quad(func_1, 0, T, args=(k,w))[0]/T,3) for k in m]# интеграл для a[k], k -номер гармоники
b=[round(2*quad(func_2, 0, T, args=(k,w))[0]/T,3) for k in m]# интеграл для b[k], k -номер гармоники
F1=[a[1]*np.cos(w*1*t)+b[1]*np.sin(w*1*t) for t in q]#функции для гармоник
F2=[a[2]*np.cos(w*2*t)+b[2]*np.sin(w*2*t) for t in q]
F3=[a[3]*np.cos(w*3*t)+b[3]*np.sin(w*3*t) for t in q]
plt.figure()
plt.title("Классический гармонический анализ функции \n при t<pi f(t)=cos(t) при 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<np.pi:
p=np.cos(t)
else:
p=-np.cos(t)
return p
def func_1(t,k,w):# функция для расчёта коэффициента a[k]
if t<np.pi:
z=np.cos(t)*np.cos(w*k*t)
else:
z=-np.cos(t)*np.cos(w*k*t)
return z
def func_2(t,k,w):#функция для расчёта коэффициента b[k]
if t<np.pi:
y=np.cos(t)*np.sin(w*k*t)
else:
y=-np.cos(t)*np.sin(w*k*t)
return y
a=[];b=[];c=4;g=[];m=np.arange(0,c,1);q=np.arange(0,2*np.pi,0.01)# подготовка списков для численного анализа
a=[round(2*quad(func_1, 0, T, args=(k,w))[0]/T,3) for k in m]# интеграл для a[k], k -номер гармоники
b=[round(2*quad(func_2, 0, T, args=(k,w))[0]/T,3) for k in m]# интеграл для b[k], k -номер гармоники
plt.figure()
plt.title("Спектральный анализ \n Спектр амплитуд-A[k]")
A=np.array([(a[k]**2+b[k]**2)**0.5 for k in m])# численные значения амплитуды гармоник
plt.plot([m[1],m[1]],[0,A[1]],label='1 гармоника')
plt.plot([m[2],m[2]],[0,A[2]],label='2 гармоника')
plt.plot([m[3],m[3]],[0,A[3]],label='3 гармоника')
plt.xlabel("Номер гармоники")
plt.ylabel("Амплитуда")
plt.legend(loc='best')
plt.grid(True)
for k in m:#вычисление численных значений фазы
if a[k]!=0:
g.append(-np.tanh(b[k]/a[k]))
else:
g.append(-np.pi/2)# фаза когда тангенс равен бесконечности
plt.figure()
plt.title("Спектральный анализ \n Спектр фаз -g(k)")
plt.plot([m[1],m[1]],[0, g[1]],label='Фаза 1 гармоники')
plt.plot([m[2],m[2]],[0, g[2]],label='Фаза 2 гармоники')
plt.plot([m[3],m[3]],[0, g[3]],label='Фаза 3 гармоники')
plt.xlabel("Номер гармоники")
plt.ylabel("Фаза")
plt.legend(loc='best')
plt.grid(True)
plt.figure()
plt.title("Спектральный синтез - FK=A[k]*cos(w*k*t+g[k])")
FK=-np.array(a[0]/2)+np.array([0*t for t in q-1])#подготовка массива длячисленного синтеза
for k in m:
FK=FK+np.array([A[k]*np.cos(w*k*t+g[k]) for t in q])# численный спектральный синтез
P=[func(t) for t in q]
plt.plot(q, P, label='f(t)')
plt.plot(q, FK, label='FK(t)')
plt.xlabel("Время t")
plt.ylabel("f(t),FK(t)")
plt.legend(loc='best')
plt.grid(True)
plt.show()
Результат
Программная реализация спектрального анализа и синтеза с использованием функций 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<np.pi:
p=np.cos(t)
else:
p=-np.cos(t)
return p
v=[f(t) for t in m]
F=np.fft.rfft(v, n=None, axis=-1) # прямое быстрое преобразование Фурье в частотную область
A=[((F[i].real)**2+(F[i].imag)**2)**0.5 for i in np.arange(0,7,1)]#модуль амплитуды
for i in np.arange(0,7,1):# определение фазы
if F[i].imag!=0:
t=(-np.tanh((F[i].real)/(F[i].imag)))
arg.append(t)
else:
arg.append(np.pi/2)
plt.figure()
plt.title("Спектральный анализ с использованием прямого БПФ ")
plt.plot(np.arange(0,7,1),arg,label='Фаза')
plt.plot(np.arange(0,7,1),A,label='Амплитуда')
plt.xlabel("Частота")
plt.ylabel("Фаза,Амплитуда")
plt.legend(loc='best')
plt.grid(True)
for i in np.arange(0,9,1):
if i<=7:
q.append(F[i])
else:
q.append(0)
h=np.fft.irfft(q, n=None, axis=-1)# обратное быстрое преобразование Фурье во временную область
plt.figure()
plt.title("Спектральный синтез с использованием обратного БПФ ")
plt.plot(m, v,label='Исходная функция')
plt.plot(m, h,label='Синтезированная функция')
plt.xlabel("Время")
plt.ylabel("Амплитуда")
plt.legend(loc='best')
plt.grid(True)
plt.show()
Фильтрация аналоговых сигналов
Под фильтрацией подразумевается выделение полезного сигнала из его смеси с мешающим сигналом — шумом. Наиболее распространенный тип фильтрации — частотная фильтрация. Если известна область частот, занимаемых полезным сигналом, достаточно выделить эту область и подавить те области, которые заняты шумом.
Программная реализация иллюстрирует технику фильтрации с применением БПФ. Сначала синтезируется исходный сигнал, представленный 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<np.pi:
p=np.cos(t)
else:
p=-np.cos(t)
return p
def FH(x):# ступенчатая функция Хэвисайда
if x>=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()
Результат
Вывод
В статье приведена попытка классификации по областям применения методов преобразования Фурье.