执行离散傅里叶变换 - FFT

Question

我正在尝试在声音处理中做一个项目，并且需要将频率放入另一个域。现在，我试图实现一个FFT，但并不顺利。我试着去了解z-transnsform，那也不好。我读了一遍，发现DFT更容易理解，尤其是算法。所以我使用示例对算法进行了编码，但我不知道或者认为输出是正确的。 （我在这里没有Matlab，也找不到任何资源来测试它），并想知道你们是否知道我是否正朝着正确的方向前进。这是我到目前为止的代码：

#include <iostream> 
#include <complex> 
#include <vector> 

using namespace std; 

const double PI = 3.141592; 

vector< complex<double> > DFT(vector< complex<double> >& theData) 
{ 
// Define the Size of the read in vector 
const int S = theData.size(); 

// Initalise new vector with size of S 
vector< complex<double> > out(S, 0); 
for(unsigned i=0; (i < S); i++) 
{ 
    out[i] = complex<double>(0.0, 0.0); 
    for(unsigned j=0; (j < S); j++) 
    { 
     out[i] += theData[j] * polar<double>(1.0, - 2 * PI * i * j/S); 
    } 
} 

return out; 
} 

int main(int argc, char *argv[]) { 

vector< complex<double> > numbers; 

numbers.push_back(102023); 
numbers.push_back(102023); 
numbers.push_back(102023); 
numbers.push_back(102023); 

vector< complex<double> > testing = DFT(numbers); 

for(unsigned i=0; (i < testing.size()); i++) 
{ 
    cout << testing[i] << endl; 

} 
} 

输入是：

102023    102023 

102023    102023 

而结果：

(408092,  0) 

(-0.0666812, -0.0666812) 

(1.30764e-07, -0.133362) 

(0.200044, -0.200043) 

任何帮助或建议将是巨大的，我不期待很多，但是，任何事情都会很棒。谢谢:)

Answer 1

您的代码看起来像是okey。 out[0]应代表输入波形的“DC”分量。在你的情况下，它比输入波形大4倍，因为你的归一化系数是1.

其他系数应该表示输入波形的幅度和相位。系数被镜像，即，out [i] == out [N-i]。你可以用下面的代码测试：

double frequency = 1; /* use other values like 2, 3, 4 etc. */ 
for (int i = 0; i < 16; i++) 
    numbers.push_back(sin((double)i/16 * frequency * 2 * PI)); 

对于frequency = 1，这给：

(6.53592e-07,0) 
(6.53592e-07,-8) 
(6.53592e-07,1.75661e-07) 
(6.53591e-07,2.70728e-07) 
(6.5359e-07,3.75466e-07) 
(6.5359e-07,4.95006e-07) 
(6.53588e-07,6.36767e-07) 
(6.53587e-07,8.12183e-07) 
(6.53584e-07,1.04006e-06) 
(6.53581e-07,1.35364e-06) 
(6.53576e-07,1.81691e-06) 
(6.53568e-07,2.56792e-06) 
(6.53553e-07,3.95615e-06) 
(6.53519e-07,7.1238e-06) 
(6.53402e-07,1.82855e-05) 
(-8.30058e-05,7.99999) 

这似乎是正确的对我说：可以忽略不计DC，幅度8月1日的谐波，其他谐波可以忽略不计的幅度。

Answer 2

@Phorce就在这里。我认为没有任何共鸣可以重新发明轮子。然而，如果你想这样做，以便你了解方法和有自己的编码喜悦，我可以提供一个FORTRAN FFT代码，我几年前开发的。当然，这不是C++，需要翻译;这应该不会太困难，应该让你学到很多东西......

下面是一个基于Radix 4的算法;这个基-4 FFT递归地将DFT分割成每四个时间样本的组的四个四分之一长度的DFT。这些更短的FFT的输出被重新用于计算许多输出，从而大大降低了总计算成本。基数-4抽取频率FFT将每第四个输出采样分组为较短长度的DFT以节省计算量。基数-4 FFT只需要基数-2 FFT的复数乘积的75％。有关更多信息，请参阅here。

!+ FILE: RADIX4.FOR 
! =================================================================== 
! Discription: Radix 4 is a descreet complex Fourier transform algorithim. It 
! is to be supplied with two real arrays, one for real parts of function 
! one for imaginary parts: It can also unscramble transformed arrays. 
! Usage: calling FASTF(XREAL,XIMAG,ISIZE,ITYPE,IFAULT); we supply the 
! following: 
! 
! XREAL - array containing real parts of transform sequence  
! XIMAG - array containing imagianry parts of transformation sequence 
! ISIZE - size of transform (ISIZE = 4*2*M) 
! ITYPE - +1 forward transform 
!   -1 reverse transform 
! IFAULT - 1 if error 
!  - 0 otherwise 
! =================================================================== 
! 
! Forward transform computes: 
!  X(k) = sum_{j=0}^{isize-1} x(j)*exp(-2ijk*pi/isize) 
! Backward computes: 
!  x(j) = (1/isize) sum_{k=0}^{isize-1} X(k)*exp(ijk*pi/isize) 
! 
! Forward followed by backwards will result in the origonal sequence! 
! 
! =================================================================== 

     SUBROUTINE FASTF(XREAL,XIMAG,ISIZE,ITYPE,IFAULT) 

     REAL*8 XREAL(*),XIMAG(*) 
     INTEGER MAX2,II,IPOW 

     PARAMETER (MAX2 = 20) 

! Check for valid transform size upto 2**(max2): 
     IFAULT = 1 
     IF(ISIZE.LT.4) THEN 
     print*,'FFT: Error: Data array < 4 - Too small!' 
     return 
     ENDIF 
     II = 4 
     IPOW = 2 

! Prepare mod 2: 
1 IF((II-ISIZE).NE.0) THEN 
     II = II*2 
     IPOW = IPOW + 1  
     IF(IPOW.GT.MAX2) THEN 
      print*,'FFT: Error: FFT1!' 
      return 
     ENDIF 
     GOTO 1 
     ENDIF 

! Check for correct type: 
     IF(IABS(ITYPE).NE.1) THEN 
     print*,'FFT: Error: Wrong type of transformation!' 
     return 
     ENDIF 

! No entry errors - continue: 
     IFAULT = 0 

! call FASTG to preform transformation: 
     CALL FASTG(XREAL,XIMAG,ISIZE,ITYPE) 

! Due to Radix 4 factorisation results are not in the same order 
! after transformation as they were when the data was submitted: 
! We now call SCRAM, to unscramble the reults: 

     CALL SCRAM(XREAL,XIMAG,ISIZE,IPOW) 

     return 

     END 

!-END: RADIX4.FOR 


! =============================================================== 
! Discription: This is the radix 4 complex descreet fast Fourier 
! transform with out unscrabling. Suitable for convolutions or other 
! applications that do not require unscrambling. Designed for use 
! with FASTF.FOR. 
! 
     SUBROUTINE FASTG(XREAL,XIMAG,N,ITYPE) 

     INTEGER N,IFACA,IFCAB,LITLA 
     INTEGER I0,I1,I2,I3 

     REAL*8 XREAL(*),XIMAG(*),BCOS,BSIN,CW1,CW2,PI 
     REAL*8 SW1,SW2,SW3,TEMPR,X1,X2,X3,XS0,XS1,XS2,XS3 
     REAL*8 Y1,Y2,Y3,YS0,YS1,YS2,YS3,Z,ZATAN,ZFLOAT,ZSIN 

     ZATAN(Z) = ATAN(Z) 
     ZFLOAT(K) = FLOAT(K) ! Real equivalent of K. 
     ZSIN(Z) = SIN(Z) 

     PI = (4.0)*ZATAN(1.0) 
     IFACA = N/4 

! Forward transform: 
     IF(ITYPE.GT.0) THEN 
     GOTO 5 
     ENDIF 

! If this is for an inverse transform - conjugate the data: 
     DO 4, K = 1,N 
     XIMAG(K) = -XIMAG(K) 
4 CONTINUE 

5 IFCAB = IFACA*4 

! Proform appropriate transformations: 
     Z = PI/ZFLOAT(IFCAB) 
     BCOS = -2.0*ZSIN(Z)**2 
     BSIN = ZSIN(2.0*Z) 
     CW1 = 1.0 
     SW1 = 0.0 

! This is the main body of radix 4 calculations: 
     DO 10, LITLA = 1,IFACA 
     DO 8, I0 = LITLA,N,IFCAB 

      I1 = I0 + IFACA 
      I2 = I1 + IFACA 
      I3 = I2 + IFACA 
      XS0 = XREAL(I0) + XREAL(I2) 
      XS1 = XREAL(I0) - XREAL(I2) 
      YS0 = XIMAG(I0) + XIMAG(I2) 
      YS1 = XIMAG(I0) - XIMAG(I2) 
      XS2 = XREAL(I1) + XREAL(I3) 
      XS3 = XREAL(I1) - XREAL(I3) 
      YS2 = XIMAG(I1) + XIMAG(I3) 
      YS3 = XIMAG(I1) - XIMAG(I3) 

      XREAL(I0) = XS0 + XS2 
      XIMAG(I0) = YS0 + YS2 

      X1 = XS1 + YS3 
      Y1 = YS1 - XS3 
      X2 = XS0 - XS2 
      Y2 = YS0 - YS2 
      X3 = XS1 - YS3 
      Y3 = YS1 + XS3 

      IF(LITLA.GT.1) THEN 
       GOTO 7 
      ENDIF 

      XREAL(I2) = X1 
      XIMAG(I2) = Y1 
      XREAL(I1) = X2 
      XIMAG(I1) = Y2 
      XREAL(I3) = X3 
      XIMAG(I3) = Y3 
      GOTO 8 

! Now IF required - we multiply by twiddle factors: 
7   XREAL(I2) = X1*CW1 + Y1*SW1 
      XIMAG(I2) = Y1*CW1 - X1*SW1 
      XREAL(I1) = X2*CW2 + Y2*SW2 
      XIMAG(I1) = Y2*CW2 - X2*SW2 
      XREAL(I3) = X3*CW3 + Y3*SW3 
      XIMAG(I3) = Y3*CW3 - X3*SW3 
8  CONTINUE 
     IF(LITLA.EQ.IFACA) THEN 
      GOTO 10 
     ENDIF 

! Calculate a new set of twiddle factors: 
     Z = CW1*BCOS - SW1*BSIN + CW1 
     SW1 = BCOS*SW1 + BSIN*CW1 + SW1 
     TEMPR = 1.5 - 0.5*(Z*Z + SW1*SW1) 
     CW1 = Z*TEMPR 
     SW1 = SW1*TEMPR   
     CW2 = CW1*CW1 - SW1*SW1 
     SW2 = 2.0*CW1*SW1 
     CW3 = CW1*CW2 - SW1*SW2 
     SW3 = CW1*SW2 + CW2*SW1 
10 CONTINUE 
     IF(IFACA.LE.1) THEN 
     GOTO 14 
     ENDIF 

! Set up tranform split for next stage: 
     IFACA = IFACA/4 
     IF(IFACA.GT.0) THEN 
     GOTO 5 
     ENDIF 

! This is the calculation of a radix two-stage: 
     DO 13, K = 1,N,2 
     TEMPR = XREAL(K) + XREAL(K + 1) 
     XREAL(K + 1) = XREAL(K) - XREAL(K + 1) 
     XREAL(K) = TEMPR 
     TEMPR = XIMAG(K) + XIMAG(K + 1) 
     XIMAG(K + 1) = XIMAG(K) - XIMAG(K + 1) 
     XIMAG(K) = TEMPR 
13 CONTINUE 
14 IF(ITYPE.GT.0) THEN 
     GOTO 17 
     ENDIF 

! For the inverse case, cojugate and scale the transform: 
     Z = 1.0/ZFLOAT(N) 
     DO 16, K = 1,N 
     XIMAG(K) = -XIMAG(K)*Z 
     XREAL(K) = XREAL(K)*Z 
16 CONTINUE 

17 return 

     END 
! ---------------------------------------------------------- 
!-END of subroutine FASTG.FOR. 
! ---------------------------------------------------------- 


!+ FILE: SCRAM.FOR 
! ========================================================== 
! Discription: Subroutine for unscrambiling FFT data: 
! ========================================================== 
     SUBROUTINE SCRAM(XREAL,XIMAG,N,IPOW) 

     INTEGER L(19),II,J1,J2,J3,J4,J5,J6,J7,J8,J9,J10,J11,J12 
     INTEGER J13,J14,J15,J16,J17,J18,J19,J20,ITOP,I 
     REAL*8 XREAL(*),XIMAG(*),TEMPR 

     EQUIVALENCE (L1,L(1)),(L2,L(2)),(L3,L(3)),(L4,L(4)) 
     EQUIVALENCE (L5,L(5)),(L6,L(6)),(L7,L(7)),(L8,L(8)) 
     EQUIVALENCE (L9,L(9)),(L10,L(10)),(L11,L(11)),(L12,L(12)) 
     EQUIVALENCE (L13,L(13)),(L14,L(14)),(L15,L(15)),(L16,L(16)) 
     EQUIVALENCE (L17,L(17)),(L18,L(18)),(L19,L(19)) 

     II = 1 
     ITOP = 2**(IPOW - 1) 
     I = 20 - IPOW 
     DO 5, K = 1,I 
     L(K) = II 
5 CONTINUE 
     L0 = II 
     I = I + 1 
     DO 6, K = I,19 
     II = II*2 
     L(K) = II 
6 CONTINUE 
     II = 0 
     DO 9, J1 = 1,L1,L0 
     DO 9, J2 = J1,L2,L1 
      DO 9, J3 = J2,L3,L2 
      DO 9, J4 = J3,L4,L3 
       DO 9, J5 = J4,L5,L4 
       DO 9, J6 = J5,L6,L5 
        DO 9, J7 = J6,L7,L6 
        DO 9, J8 = J7,L8,L7 
         DO 9, J9 = J8,L9,L8 
         DO 9, J10 = J9,L10,L9 
          DO 9, J11 = J10,L11,L10 
          DO 9, J12 = J11,L12,L11 
           DO 9, J13 = J12,L13,L12 
           DO 9, J14 = J13,L14,L13 
            DO 9, J15 = J14,L15,L14 
            DO 9, J16 = J15,L16,L15 
             DO 9, J17 = J16,L17,L16 
             DO 9, J18 = J17,L18,L17 
              DO 9, J19 = J18,L19,L18 
              J20 = J19 
              DO 9, I = 1,2 
               II = II +1 
               IF(II.GE.J20) THEN 
                GOTO 8 
               ENDIF 
! J20 is the bit reverse of II! 
! Pairwise exchange: 
               TEMPR = XREAL(II) 
               XREAL(II) = XREAL(J20) 
               XREAL(J20) = TEMPR 
               TEMPR = XIMAG(II) 
               XIMAG(II) = XIMAG(J20) 
               XIMAG(J20) = TEMPR 
8            J20 = J20 + ITOP 
9 CONTINUE 

     return 

     END 
! ------------------------------------------------------------------- 
!-END: 
! ------------------------------------------------------------------- 

经过这个并理解它需要时间！我使用多年前发现的CalTech论文写了这篇文章，我无法回想起恐惧的参考文献。祝你好运。

我希望这会有所帮助。

Answer 3

您的代码有效。 我会给更多数字的PI（3.1415926535898）。 此外，您必须将DFT总和的输出与D的大小相除。

由于测试中的输入序列是恒定的，因此DFT输出应该只有一个非零系数。 事实上，所有的输出系数都相对于第一个系数非常小。

但是对于较大的输入长度，这不是实现DFT的有效方式。 如果计时是一个问题，请查看Fast Fourrier Transform以获得更快速的计算DFT的方法。

Answer 4

您的代码看起来很适合我。我不确定你对输出的期望是什么，但是，考虑到你的输入是一个常数值，一个常数的DFT是一个在仓0中的DC项，在其余的仓中为零（或者你有一个相近的等价物） 。

您可能会尝试使用包含某种类型的波形（如正弦波或方波）的较长序列来测试您的代码。但是，一般来说，你应该考虑在生产代码中使用类似fftw的东西。很长一段时间以来，它一直被很多人所忽视和高度优化。对于特殊情况（例如，2的幂的长度），FFT是优化的DFT。

Answer 5

MoonKnight已经在Fortran中提供了一个基数-4抽取频率 Cooley-Tukey方案。我在下面提供了一个radix-2 Decimation In Frequency Cooley-Tukey scheme in Matlab。

该代码是一个迭代之一，考虑在以下图中的方案：

递归方法也是可能的。如您所见，该实现还会计算执行的乘法和加法运算的次数，并将其与How many FLOPS for FFT?中报告的理论计算结果进行比较。该代码明显比由Matlab利用的高度优化的FFTW慢得多。

还要注意，旋转因子omegaa^((2^(p - 1) * n))可以脱机计算，然后从查找表中恢复，但是在下面的代码中跳过了这一点。

对于matlab实现迭代基2抽取在时间库里 - 杜克方案，请参阅Implementing a Fast Fourier Transform for Option Pricing。

% --- Radix-2 Decimation In Frequency - Iterative approach 

clear all 
close all 
clc 

N = 32; 

x = randn(1, N); 
xoriginal = x; 
xhat = zeros(1, N); 

numStages = log2(N); 

omegaa = exp(-1i * 2 * pi/N); 

mulCount = 0; 
sumCount = 0; 
tic 
M = N/2; 
for p = 1 : numStages; 
    for index = 0 : (N/(2^(p - 1))) : (N - 1); 
     for n = 0 : M - 1;   
      a = x(n + index + 1) + x(n + index + M + 1); 
      b = (x(n + index + 1) - x(n + index + M + 1)) .* omegaa^((2^(p - 1) * n)); 
      x(n + 1 + index) = a; 
      x(n + M + 1 + index) = b; 
      mulCount = mulCount + 4; 
      sumCount = sumCount + 6; 
     end; 
    end; 
    M = M/2; 
end 
xhat = bitrevorder(x); 
timeCooleyTukey = toc; 

tic 
xhatcheck = fft(xoriginal); 
timeFFTW = toc; 

rms = 100 * sqrt(sum(sum(abs(xhat - xhatcheck).^2))/sum(sum(abs(xhat).^2))); 

fprintf('Time Cooley-Tukey = %f; \t Time FFTW = %f\n\n', timeCooleyTukey, timeFFTW); 
fprintf('Theoretical multiplications count \t = %i; \t Actual multiplications count \t = %i\n', ... 
     2 * N * log2(N), mulCount); 
fprintf('Theoretical additions count \t\t = %i; \t Actual additions count \t\t = %i\n\n', ... 
     3 * N * log2(N), sumCount); 
fprintf('Root mean square with FFTW implementation = %.10e\n', rms);