从离散傅里叶变换建模傅里叶级数进行外推

Question

我正在尝试将 python numpy/scipy 的 fft、rfft 和 dct 转换回 sine/cosine 波的总和以重建原始数据集。我想这样做是因为我希望能够用 more/less 个采样点重建原始数据集（我相信它可能已经被 scipy.signal.resample 覆盖）并且主要是因为我想扩展 sine/cosine 未来的系列与在某些系列中如何使用线性回归来给出未来值的一般概念没有什么不同。我知道这在技术上是不正确的，因为 fft 假设离散样本在所有时间点重复，dct 假设数据是 "mirrored"，但我认为它可能具有某种短期预测值。

我试着按照这里写的内容作为 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()

这段代码应该做的是：

将所有输入数据数组转换为复数（因为 dct 不输出复数数据类型）。
归一化傅立叶系数，因为 fft() 和相关变换似乎不会除以数据集中的元素数。
填充代表个体的 lambda 函数数组每个傅里叶频率的贡献（我假设它们是傅里叶系数的序数位置）
对每个正弦曲线的个体贡献进行累计求和在要采样以重建的新点处采用的 lambda 函数一个系列

在此特定代码中，我试图查看我的重组是否等于原始 series/scipy 的逆分解，只是为了确保我做对了。我认为代码工作正常，但它用于 sine/cosine 重建的基础公式是错误的。这是此特定代码的输出：

绿色是我重建的值，Orange/Blue是原始值。显然，我的算法没有正确地重新制作系列。使用振幅和相位将正弦项和余弦项组合成一个余弦项，正如其他网站所推荐的那样，给出了不同但仍然不正确的结果，很可能是由于减去上面推荐的正弦项资源。有谁知道我的公式或代码是怎么错的？我认为它要么在 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 在这种方式是周期性的。即y[N]==y[0]、y[N+1]==y[1]等

从离散傅里叶变换建模傅里叶级数进行外推

Modeling a Fourier Series from Discrete Fourier Transform for Extrapolation

python

signal-processing

fft

dct

scipy