2

我正在尝试做一个声音处理项目，需要将频率放入另一个域。现在，我尝试实现 FFT，但效果并不好。我试图理解z-transform，但也不太顺利。我阅读并发现 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)

任何帮助或建议都会很棒，我并不期待很多，但是，任何事情都会很棒。谢谢 ：）

4

6 回答 6

6

@Phorce 就在这里。我认为没有任何理由重新发明轮子。但是，如果您想这样做以便理解该方法并享受自己编写代码的乐趣，我可以提供我几年前开发的 FORTRAN FFT 代码。当然这不是 C++，需要翻译；这应该不会太难，应该能让你学到很多东西......

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

!+ 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:
! -------------------------------------------------------------------

经历这些并理解它需要时间！我使用几年前发现的加州理工学院论文写了这篇文章，我恐怕不记得参考了。祝你好运。

我希望这有帮助。

于 2012-09-03T16:50:22.793 回答
2

您的代码有效。我会为 PI 提供更多数字（ 3.1415926535898 ）。此外，您必须将 DFT 求和的输出除以 S，即 DFT 大小。

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

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

于 2012-09-03T17:05:51.483 回答
2

你的代码对我来说很合适。我不确定您对输出的期望是什么，但是，鉴于您的输入是一个常数值，常数的 DFT 是 bin 0 中的 DC 项，其余 bin 中为零（或您拥有的接近等价物） .

您可以尝试使用包含某种波形（如正弦波或方波）的较长序列来测试您的代码。但是，一般来说，您应该考虑在生产代码中使用类似 fftw 的东西。长期以来，它被许多人绞尽脑汁并高度优化。FFT 是针对特殊情况优化的 DFT（例如，长度为 2 的幂）。

于 2012-09-03T17:13:08.293 回答
1

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

其他系数应代表输入波形的幅度和相位。系数是镜像的，即 out[i] == out[Ni]。您可以使用以下代码对此进行测试：

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)

这对我来说似乎是正确的：直流可忽略不计，第一次谐波的幅度为 8，其他谐波的幅度可忽略不计。

于 2012-09-03T17:02:45.133 回答
0

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

代码是迭代的，考虑下图中的方案：

在此处输入图像描述

递归方法也是可能的。

正如您将看到的，该实现还计算了执行的乘法和加法的数量，并将其与FFT 的多少 FLOPS 中报告的理论计算进行比较？.

该代码显然比 Matlab 使用的高度优化的 FFTW 慢得多。

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

有关迭代radix-2 时间抽取Cooley-Tukey 方案的 Matlab 实现，请参阅实现期权定价的快速傅里叶变换。

% --- 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);
于 2017-02-18T07:44:24.693 回答
0

您的代码对于获取 DFT 是正确的。

您正在测试的函数(sin ((double) i / points * frequency * 2)对应于幅度为 1、频率为 1 和采样频率 Fs = 所取点数的同步曲线。

使用获得的数据进行操作，我们有：

在此处输入图像描述

如您所见，DFT 系数关于位置系数 N / 2 是对称的，因此只有前 N / 2 提供信息。通过实部和虚部的模获得的幅度必须除以N并乘以2才能重建。系数的频率将是 Fs / N 乘以系数数的倍数。

如果我们引入两个正弦曲线，一个是频率 2 和幅度 1.3，另一个是频率 3 和幅度 1.7。

for (int i = 0; i < 16; i++)
{
    numbers.push_back(1.3 *sin((double)i / 16 * frequency1 * 2 * PI)+ 1.7 *
        sin((double)i / 16 * frequency2 * 2 * PI));
}

得到的数据是：

在此处输入图像描述

祝你好运。

于 2018-05-13T09:05:10.297 回答