Данные ряда Фурье соответствуют NumPy: FFT против кодирования - PullRequest
Новая
       16

Данные ряда Фурье соответствуют NumPy: FFT против кодирования

0 голосов
/ 11 февраля 2019

Предположим, у меня есть некоторые данные, y, к которым я хотел бы привести ряд Фурье.На этой записи было опубликовано решение Mermoz с использованием сложного формата ряда и «вычисления коэффициента с помощью суммы Римана».На этом другом посте серия получена через БПФ, и записан пример.

Я пытался реализовать оба подхода (изображение и код ниже - обратите внимание, каждый раз, когда код запускается, разныеданные будут генерироваться из-за использования numpy.random.normal ), но мне интересно, почему я получаю разные результаты - подход Римана кажется «ошибочно смещенным», в то время как подход БПФ кажется «сжатым».Я также не уверен в моем определении периода "тау" для ряда.Я ценю внимание.

Я использую Spyder с Python 3.7.1 в Windows 7

Пример 

import matplotlib.pyplot as plt
import numpy as np

# Assume x (independent variable) and y are the data.
# Arbitrary numerical values for question purposes:
start = 0
stop = 4
mean = 1
sigma = 2
N = 200
terms = 30 # number of terms for the Fourier series

x = np.linspace(start,stop,N,endpoint=True) 
y = np.random.normal(mean, sigma, len(x))

# Fourier series
tau = (max(x)-min(x)) # assume that signal length = 1 period (tau)

# From ref 1
def cn(n):
    c = y*np.exp(-1j*2*n*np.pi*x/tau)
    return c.sum()/c.size
def f(x, Nh):
    f = np.array([2*cn(i)*np.exp(1j*2*i*np.pi*x/tau) for i in range(1,Nh+1)])
    return f.sum()
y_Fourier_1 = np.array([f(t,terms).real for t in x])

# From ref 2
Y = np.fft.fft(y)
np.put(Y, range(terms+1, len(y)), 0.0) # zero-ing coefficients above "terms"
y_Fourier_2 = np.fft.ifft(Y)

# Visualization
f, ax = plt.subplots()
ax.plot(x,y, color='lightblue', label = 'artificial data')
ax.plot(x, y_Fourier_1, label = ("'Riemann' series fit (%d terms)" % terms))
ax.plot(x,y_Fourier_2, label = ("'FFT' series fit (%d terms)" % terms))
ax.grid(True, color='dimgray', linestyle='--', linewidth=0.5)
ax.set_axisbelow(True)
ax.set_ylabel('y')
ax.set_xlabel('x')
ax.legend()

1 Ответ

0 голосов
/ 12 февраля 2019

Выполнение двух небольших модификаций достаточно, чтобы суммы были почти аналогичны выводу np.fft. Библиотека FFTW действительно вычисляет эти суммы .

1) Среднее значение сигнала c[0] должно учитываться: 

f = np.array([2*cn(i)*np.exp(1j*2*i*np.pi*x/tau) for i in range(0,Nh+1)]) # here : 0, not 1

2) Выход должен быть масштабирован.

y_Fourier_1=y_Fourier_1*0.5

enter image description here Выход выглядит «сжатым», поскольку высокочастотные компоненты были отфильтрованы.Действительно, высокочастотные колебания на входе были очищены, а выход выглядит как скользящее среднее.

Здесь tau фактически определяется как stop-start: оно соответствует длине кадра.Это ожидаемый период сигнала.

Если кадр не соответствует периоду сигнала, вы можете угадать его период, свернув сигнал с самим собой и найдя первый максимум.См. Поиск периода сигнала вне БПФ Тем не менее, вряд ли он будет правильно работать с набором данных, сгенерированным numpy.random.normal: это аддитивный белый гауссовский шум .Поскольку он имеет постоянную спектральную плотность мощности, его вряд ли можно считать периодическим!

...