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

Question

使用 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,end=100), col="red",   type="l",lty=1,lwd=1, xaxt="n");
   axis(1, at=seq(round(t[1],1),round(t[length(t)],1), 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 算法的问题的了解都会有所帮助。我不想只是黑箱我对函数的使用 - 我想在继续更复杂的函数之前更深入地理解这些事情。

谢谢， 克里斯蒂安

Answer 1

主要问题来自信号索引。首先要获得 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=0

and

k=0` 的参数，您需要调整索引：

kk=seq(1,N, 1);  nn=seq(1,N, 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,N, 1);  nn=seq(1,N, 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,end=100), col="red",   type="l",lty=2,lwd=1, xaxt="n");
axis(1, at=seq(round(t[1],1),round(t[length(t)],1), by=0.1), las=2);
max(abs(sig.ts - sig.freqd_inv)) / max(sig.ts)

这应该会产生一个围绕

1.814886e-13

的指标，这可能更符合您的预期。相应的图还应该显示原始信号和往返信号重叠：