Моделирование ряда Фурье из дискретного преобразования Фурье для экстраполяции

Question

Я пытаюсь обратить преобразования Python Numpy / Scipy's fft, rfft и dct обратно в сумму синусоидальных / косинусных волн, чтобы восстановить исходный набор данных. Я хочу сделать это, потому что я хочу иметь возможность восстановить исходный набор данных с большим / меньшим количеством точек выборки (которые, я считаю, уже могут быть покрыты scipy.signal.resample) и главным образом потому, что я хочу расширить ряд синус / косинус в будущее мало чем отличается от того, как линейная регрессия может использоваться в некоторых сериях, чтобы дать общее представление о будущих ценностях. Я знаю, что это технически неверно, так как fft предполагает, что дискретная выборка повторяется во всех временных точках, а dct предполагает, что данные «зеркально отражены», но я думаю, что они могут иметь какое-то краткосрочное прогностическое значение.

Я пытался следовать тому, что здесь было написано, как руководство к алгоритмам декомпозиции Numpy: http://snowball.millersville.edu/~adecaria/ESCI386P/esci386-lesson17-Fourier-Transforms.pdf

Вот мой код:

import numpy as np
from scipy.fftpack import fft,ifft,dct,idct,rfft,irfft
import matplotlib.pyplot as plt

def reconstructSeries(transformedVals,newxvals):
    transformedVals=transformedVals.astype('complex128')

    transformedVals=transformedVals/len(transformedVals) #for some reason, numpy does not normalize the values it has, so I have to do it here.


    reconstructedVals=np.zeros(len(newxvals))
    series=[]
    # perhaps [:len(transformedVals)//2] ?
    for frequency,val in enumerate(transformedVals):    
        #the position of the coefficient is the frequency (in radians)

        #amplitude=np.sqrt(np.real(val)**2+np.imag(val)**2)
        #phase=np.arctan(np.imag(val)/np.real(val))
        series.append(lambda x: np.real(val)*np.cos(frequency*newxvals)-np.imag(val)*np.sin(frequency*newxvals))
        #series.append(lambda x: amplitude*np.cos(2*np.pi*frequency*newxvals+phase)) #this is in radians to accomidate phase and the default cosine function
        reconstructedVals=reconstructedVals+np.array(series[frequency](newxvals))

    return reconstructedVals,series

#y=np.arange(250)
y=np.cos(np.arange(250)+5)

yf = fft(y) #this can be rfft or dct as well
myyvalues,sinosoidseries=reconstructSeries(yf,np.arange(250))

plt.plot(ifft(yf));plt.plot(y);plt.plot(myyvalues);plt.show()

Что должен делать этот код:

1. Преобразовать все входные массивы данных в сложный (поскольку dct не выводит сложный тип данных).
2. Нормализуйте коэффициенты Фурье, так как fft () и связанные преобразования не делятся на количество элементов в наборе данных.
3. Заполнить массив лямбда-функций, представляющих человека вклад каждой частоты Фурье (я предполагаю, что они являются порядковыми позициями коэффициентов Фурье)
4. Совокупная сумма отдельных вкладов каждого синусоидального лямбда-функция, взятая в новых точках для выборки для реконструкции серия

В этом конкретном коде я пытаюсь увидеть, равна ли моя перекомпоновка обратной декомпозиции оригинальной серии / scipy, просто чтобы убедиться, что я все делаю правильно. Я думаю, что код работает нормально, но основные формулы, которые он использует для восстановления синуса / косинуса, неверны. Вот вывод для этого конкретного кода:

Зеленый - мои восстановленные значения, а оранжевый / синий - исходные значения. Ясно, что мой алгоритм не корректирует серию должным образом. Использование амплитуды и фазы для объединения терминов синуса и косинуса в один термин косинуса, как было рекомендовано на других сайтах, дает другой, но все же неверный результат, скорее всего из-за вычитания термина синуса, рекомендованного в приведенном выше источник. Кто-нибудь знает, как моя формула или код неверен? Я думаю, что это либо в части cos () - sin (), либо что-то вроде частоты, не умноженной на константу.

* Примечание: я знаю, что этот вопрос похож на: ряд Фурье из дискретного преобразования Фурье но я не думаю, что ответ на этот вопрос мне подходит здесь.

Answer 1

Ошибка, которую я вижу в коде, заключается в сложном умножении: вы умножаете действительный компонент выборки в частотной области на cos, а мнимый компонент на sin. Это не то, как работает сложное умножение. Вам нужно умножить комплексное значение выборки на комплексное значение cos + i sin. Комплексные числа a + ib и c + id при умножении дают ac-bd + iad + ibc, а не ac + bd.

---

Редактировать: Как заполнить частотную область нулями для интерполяции

Функция SciPy ifft имеет параметр n, который можно использовать для заполнения массива нулями перед преобразованием. Не используйте этот параметр. Он добавляет нули к концу сигнала, нарушая симметрию сигнала, и, следовательно, обычно дает нереальный результат.

DFT (то, что fft вычисляет) имеет частоты k = 0 ... N-1. Но k является периодическим, что означает, что k=N-1 совпадает с k=-1. Нам нужно сохранить (комплексно сопряженную) симметрию вокруг k=0, которая существует в области Фурье для реальных значений сигналов во временной области, что означает, что значения сигнала в частотной области для k=1 и k=-1 должны поддерживаться эта симметрия (а также для k=2 и k=-2 и т. д.).

При заполнении нулями мы увеличиваем это значение N, поэтому мы также меняем местоположение k=-1 в массиве (как и в k=N-1, увеличение N означает, что это местоположение перемещается).

Таким образом, отступы должны добавлять нули прямо в середине массива, чтобы исходные значения в начале и конце массива сохранялись. Середина массива находится между (N+1)//2-1 и (N+1)//2:

N = 250
y = np.cos(np.arange(N)+5)
yf = fft(y)
yf = np.concatenate((yf[:(N+1)//2], np.zeros(N), yf[(N+1)//2:]))
y2 = ifft(yf)
plt.subplot(2,1,1)
plt.plot(y,'.-')
plt.subplot(2,1,2)
plt.plot(y2,'.-')
plt.show()

Обратите внимание, что сигнал во временной области остается неизменным, но имеет удвоенное количество выборок.

Обратите внимание также, как это не экстраполирует: если вы расширяете синусоидальные и косинусные волны, составляющие y, вы будете восстанавливать значения с начала сигнала, так как y, реконструируемый таким образом, является периодическим. То есть y[N]==y[0], y[N+1]==y[1] и т. Д.