c++ - 離散フーリエ変換の実装-FFT

Question

私はサウンド処理のプロジェクトを行おうとしていますが、周波数を別のドメインに入れる必要があります。今、私はFFTを実装しようとしましたが、うまくいきませんでした。私は-transformを理解しようとしましたがz、それもうまくいきませんでした。私が読んだところ、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

@Phorceはここにあります。車輪の再発明をする理由はないと思います。ただし、方法論を理解し、自分でコーディングする喜びを味わうためにこれを実行したい場合は、数年前に開発したFORTRANFFTコードを提供できます。もちろん、これはC ++ではなく、翻訳が必要になります。これはそれほど難しいことではなく、そうすることで多くを学ぶことができるはずです...

以下は、Radix4ベースのアルゴリズムです。この基数4FFTは、DFTを4回ごとの時間サンプルのグループの4つの1/4長DFTに再帰的に分割します。これらの短いFFTの出力は、多くの出力を計算するために再利用されるため、総計算コストが大幅に削減されます。基数4の周波数デシメーションFFTは、4つおきの出力サンプルをより短い長さのDFTにグループ化して、計算を節約します。基数4のFFTは、基数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:
! -------------------------------------------------------------------

これを通り抜けて理解するには時間がかかります！私は何年も前に見つけたCalTechの論文を使ってこれを書きました、私は恐れている参照を思い出すことができません。幸運を。

これがお役に立てば幸いです。

Answer 2

コードは機能します。PI（3.1415926535898）の桁数を増やします。また、DFTの合計の出力をS（DFTサイズ）で割る必要があります。

テストの入力系列は一定であるため、DFT出力にはゼロ以外の係数が1つだけ含まれている必要があります。そして実際、すべての出力係数は最初の係数に比べて非常に小さいです。

ただし、入力長が長い場合、これはDFTを実装する効率的な方法ではありません。タイミングが懸念される場合は、高速フーリエ変換を調べて、DFTを計算するためのより高速な方法を探してください。

Answer 3

あなたのコードは私には正しく見えます。出力に何を期待していたかはわかりませんが、入力が定数値であるとすると、定数のDFTは、ビン0ではDC項であり、残りのビン（またはそれに近いもの）ではゼロです。 。

正弦波や方形波などのある種の波形を含む、より長いシーケンスでコードをテストしてみてください。ただし、一般的には、本番コードでfftwのようなものを使用することを検討する必要があります。それは長い間多くの人々によって絞り出され、高度に最適化されています。FFTは、特殊なケース（たとえば、2の累乗の長さ）用に最適化されたDFTです。

Answer 4

あなたのコードは大丈夫に見えます。out[0]入力波形の「DC」成分を表す必要があります。あなたの場合、正規化係数が1であるため、入力波形の4倍の大きさになります。

他の係数は、入力波形の振幅と位相を表す必要があります。係数はミラーリングされます。つまり、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)

これは私には正しいように思われます：無視できるDC、第1高調波の振幅8、他の高調波の振幅は無視できます。

Answer 5

MoonKnightは、Fortranで基数4の周波数Cooley-Tukeyスキームをすでに提供しています。以下では、Matlabで基数2のDecimation InFrequencyCooley -Tukeyスキームを提供しています。

このコードは反復的なものであり、次の図のスキームを考慮しています。

再帰的なアプローチも可能です。

ご覧のとおり、実装では、実行された乗算と加算の数も計算され、「FFTのFLOPSの数」で報告されている理論計算と比較されます。。

このコードは、Matlabが利用する高度に最適化されたFFTWよりも明らかにはるかに低速です。

また、ツイドルファクターomegaa^((2^(p - 1) * n))はオフラインで計算してルックアップテーブルから復元することもできますが、以下のコードではこの点をスキップしています。

反復基数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);

Answer 6

あなたのコードはDFTを取得するために正しいです。

テストしている関数は(sin ((double) i / points * frequency * 2)、振幅1、周波数1、およびサンプリング周波数Fs=取得したポイントの数のシノイドに対応します。

得られたデータを使用して操作します。

ご覧のとおり、DFT係数は位置係数N / 2に関して対称であるため、最初のN/2のみが情報を提供します。実数部と虚数部のモジュールによって得られた振幅は、それを再構築するためにNで除算し、2を掛ける必要があります。係数の周波数は、係数数によるFs/Nの倍数になります。

2つの正弦波を導入すると、1つは周波数2と振幅1.3、もう1つは周波数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));
} 

得られたデータは次のとおりです。

幸運を。