基于离散傅立叶变换的傅里叶级数的外推建模

Question

我试图将python numpy/scipy的fft、rfft和dct变换转换回正弦波/余弦波的总和，以重建原始数据集。我想这样做是因为我想用更多/更少的采样点重建原始数据集（我相信可能已经包含了scipy.signal.重采样)主要是因为我想把正弦/余弦级数扩展到未来，就像线性回归在某些序列中的使用一样，给出未来价值。我知道这在技术上是不正确的，因为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. 将所有输入数据数组转换为complex（因为dct不输出复杂的数据类型）。在
2. 规范化傅立叶系数，因为fft（）和相关变换似乎没有除以数据集中的元素数。在
3. 填充表示个体的lambda函数的数组 每个傅里叶频率的贡献（我假设它们是傅里叶系数的序数位置）
4. 将每个正弦曲线的单独贡献累加起来 在新采样点取lambda函数进行重构 系列

在这段特定的代码中，我试图查看我的重组是否等于原始series/scipy的逆分解，以确保我做得正确。我认为代码运行良好，但它用于正弦/余弦重建的基本公式是错误的。以下是此特定代码的输出：

绿色是我重建的值，橙色/蓝色是原始值。显然，我的算法没有正确地重新生成序列。使用振幅和相位将正弦和余弦项合并为一个余弦项，正如在其他站点上所建议的那样，给出了不同但仍然不正确的结果，很可能是由于上述来源中建议的正弦项被减去。有人知道我的公式或代码是怎么错的吗？我想它要么在cos（）-sin（）部分，要么是频率不乘以常数。在

*注：我知道这个问题有点像： Fourier Series from Discrete Fourier Transform 但我不认为这个问题的答案对我有效。在

Answer 1

我在代码中看到的错误是复数乘法：频域样本的实分量乘以cos，虚分量乘以sin。这不是复数乘法的工作原理。您需要将复采样值乘以复值cos+isin。复数a+ib和c+id相乘得到ac bd+iad+ibc，而不是ac+bd

---

编辑：如何用零点填充频域进行插值

The SciPy ^{} function有一个参数n，可用于在转换之前用零填充数组。不要使用此参数。它在信号的末尾加零，破坏了信号的对称性，因此通常会产生一个非实数的结果。在

DFT（由fft计算的内容）的频率是k=0…N-1。但是k是周期性的，这意味着{}与{}是相同的。对于实值时域信号，我们需要保持存在于Fourier域中的k=0周围的（复共轭）对称性，这意味着k=1和{}的频域信号值必须保持这种对称性（以及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]，等等