R：逆 fft() 以确认我的手动 DFT 算法不准确？

问题描述

使用 R，在评估我自己手动实施 DFT 的某些准确性指标之前，我想通过执行以下操作对 stats::fft() 的执行情况进行完整性检查：

sig.ts = ts( sin(2*pi*freq1*t) + sin(2*pi*freq2*t) );
sig.rt = fft(fft(sig.ts)/N,inverse="true");
#the two plots so perfectly align that you can't see them both
max(abs(sig.ts - sig.rt)) / max(sig.ts);
#arbitrary crude accuracy metric=1.230e-15 - EXCELLENT!

但我想自己为 DFT 编写代码，以确保我理解它，然后将其反转以希望它是相同的：

##The following is the slow DFT for Now,not the FFT...
sR = 102.4;  #the number of Hz at which we sample
freq1=3; freq2=12;  #frequency(ies) of the wave
t = seq(1/sR,10,1/sR);
sig.ts = ts( sin(2*pi*freq1*t) + sin(2*pi*freq2*t) );
N=length(t);  kk=seq(0,N/2-1,1);  nn=seq(0,N-1,1);
for(k in kk){ 
  sig.freqd[k]=0;
  for(n in nn){
    sig.freqd[k] = sig.freqd[k] + sig.ts[n+1]*exp(-j*2*pi*n*k/N);  } }
sig.freqd = (1/N)*sig.freqd; #for normalization

#Checking the "accuracy" of my manual implementation of DFT...
sig.freqd_inv=Re(fft(sig.freqd,inverse="true"));
plot(t[1:100],window(sig.ts,end=100),col="black",type="l",lty=1,lwd=1,xaxt="n");  
lines(t[1:100],window(sig.freqd_inv,col="red",xaxt="n");
   axis(1,at=seq(round(t[1],1),round(t[length(t)],by=0.1),las=2);
max(abs(sig.ts[1:(N/2-1)] - sig.freqd_inv)) / max(sig.ts[1:(N/2-1)]);  #the metric here =1.482 unfortunately

即使没有指标，该图也很明显地表明这里有些问题 - 幅度较低，可能异相，而且锯齿状。在我所有的自学中，我会说我对这一切对向量长度的敏感程度有点困惑......以及如何确保在绘图时考虑虚部的相位信息。

最重要的是，任何对我的 DFT 算法问题的见解都会有所帮助。我不想只是黑箱我对函数的使用 - 我想在继续使用更复杂的函数之前更深入地了解这些东西。

谢谢，克里斯蒂安

解决方法

主要问题来自信号索引。首先要获得 R 的 fft(...,inverse = TRUE) 可用的完整变换，您需要计算所有 N 系数（即使 N/2-1 以上的系数可以通过对称获得）。然后你应该意识到 R 中的数组索引是基于 1 的。因此，在索引 sig.freqd[k] 时，索引 k 应该从 1 而不是 0 开始。由于 exp(-1i*2*pi*n*k/N) should start with n=0andk=0` 的参数，您将需要调整指数：

kk=seq(1,N,1);  nn=seq(1,1);
for(k in kk){ 
  sig.freqd[k]=0i;
  for(n in nn){
    sig.freqd[k] = sig.freqd[k] + sig.ts[n]*exp(-1i*2*pi*(n-1)*(k-1)/N);
  }
}

我还更改了您使用 j 来表示虚数 1i 的用法，因为这是 R 识别的常用符号（并且 R 在按原样尝试您发布的示例时抱怨它） .如果您定义了不应影响结果的 j=1i。

还要注意 R 的 fft 是非标准化的。因此，为了获得相同的正向变换结果，您的 DFT 实现不应包括 1/N 规范化。另一方面，您需要添加这个因子作为最后一步，以获得与原始信号匹配的全圆正向+反向变换。

通过这些更改，您应该拥有以下代码：

##The following is the slow DFT for now,not the FFT...
sR = 102.4;  #the number of Hz at which we sample
freq1=3; freq2=12;  #frequency(ies) of the wave
t = seq(1/sR,10,1/sR);
sig.ts = ts( sin(2*pi*freq1*t) + sin(2*pi*freq2*t) );
N=length(t);  kk=seq(1,1);
for(k in kk){ 
  sig.freqd[k]=0i;
  for(n in nn){
    sig.freqd[k] = sig.freqd[k] + sig.ts[n]*exp(-1i*2*pi*(n-1)*(k-1)/N);
  }
}

#Checking the "accuracy" of my manual implementation of DFT...
sig.freqd_inv=(1/N)*Re(fft(sig.freqd,inverse="true"));
plot(t[1:100],window(sig.ts,end=100),col="black",type="l",lty=1,lwd=2,xaxt="n");  
lines(t[1:100],window(sig.freqd_inv,col="red",lty=2,lwd=1,xaxt="n");
axis(1,at=seq(round(t[1],1),round(t[length(t)],by=0.1),las=2);
max(abs(sig.ts - sig.freqd_inv)) / max(sig.ts)

这应该会产生一个围绕 1.814886e-13 的指标，这可能更符合您的预期。相应的图还应显示原始信号和往返信号重叠：

dft fft fft r r spectral