python - 通过离散傅里叶变换对傅里叶级数进行建模以进行外推

Question

我正在尝试将 python numpy/scipy 的 fft、rfft 和 dct 转换回正弦/余弦波之和，以重建原始数据集。我想这样做是因为我希望能够用更多/更少的采样点重建原始数据集(我相信可能已经被 scipy.signal.resample 覆盖)，主要是因为我想将正弦/余弦系列扩展到 future 与在某些系列中使用线性回归来给出 future 值的一般概念没有什么不同。我知道这在技术上是不正确的，因为 fft 假设离散样本在所有时间点上重复，而 dct 假设数据是“镜像的”，但我认为它可能具有某种短期预测值(value)。

我尝试遵循此处编写的内容作为 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. 填充代表个体的 lambda 函数数组每个傅里叶频率的贡献(我假设它们是傅里叶系数的顺序位置)
4. 将每个正弦曲线的单独贡献相加在要采样以重建的新点处采用的 lambda 函数一个系列

在这个特定的代码中，我试图看看我的重组是否等于原始系列/scipy 的逆分解，只是为了确保我做得正确。我认为代码工作正常，但它用于正弦/余弦重建的基础公式是错误的。以下是此特定代码的输出:

绿色是我重建的值，橙色/蓝色是原始值。显然，我的算法没有正确地重新制作该系列。正如其他网站上所建议的那样，使用幅度和相位将正弦和余弦项组合成单个余弦项，会给出不同但仍然不正确的结果，很可能是由于减去了上面建议的正弦项来源。有谁知道我的公式或代码是错误的吗？我认为它要么是在 cos()-sin() 部分，要么是频率没有乘以常数之类的东西。

*注意:我知道这个问题有点像: Fourier Series from Discrete Fourier Transform但我认为这个问题的答案对我来说不起作用。

Answer 1

我在代码中看到的错误在于复数乘法:将频域样本的实部与 cos 相乘，将虚部与 sin 相乘。这不是复数乘法的工作原理。您需要将复数样本值与复数值 cos + i sin 相乘。复数 a+ib 和 c+id 相乘时产生 ac-bd+iad+ibc，而不是 ac+bd。

<小时/>

编辑:如何用零填充频域以进行插值

The SciPy ifft function有一个参数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以这种方式重构的code>是周期性的。即y[N]==y[0]、y[N+1]==y[1]等