Implementation of the Discrete Fourier Transform - FFT

@Phorce is right here. I don't think there is any reson to reinvent the wheel. However, if you want to do this so that you understand the methodology and to have the joy of coding it yourself I can provide a FORTRAN FFT code that I developed some years ago. Of course this is not C++ and will require a translation; this should not be too difficult and should enable you to learn a lot in doing so...

Below is a Radix 4 based algorithm; this radix-4 FFT recursively partitions a DFT into four quarter-length DFTs of groups of every fourth time sample. The outputs of these shorter FFTs are reused to compute many outputs, thus greatly reducing the total computational cost. The radix-4 decimation-in-frequency FFT groups every fourth output sample into shorter-length DFTs to save computations. The radix-4 FFTs require only 75% as many complex multiplies as the radix-2 FFTs. See here for more information.

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

Going through this and understanding it will take time! I wrote this using a CalTech paper I found years ago, I cannot recall the reference I am afraid. Good luck.

I hope this helps.

Your code works. I would give more digits for PI ( 3.1415926535898 ). Also, you have to devide the output of the DFT summation by S, the DFT size.

Since the input series in your test is constant, the DFT output should have only one non-zero coefficient. And indeed all the output coefficients are very small relative to the first one.

But for a large input length, this is not an efficient way of implementing the DFT. If timing is a concern, look into the Fast Fourrier Transform for faster methods to calculate the DFT.

Your code looks right to me. I'm not sure what you were expecting for output but, given that your input is a constant value, the DFT of a constant is a DC term in bin 0 and zeroes in the remaining bins (or a close equivalent, which you have).

You might try testing you code with a longer sequence containing some type of waveform like a sine wave or a square wave. In general, however, you should consider using something like fftw in production code. Its been wrung out and highly optimized by many people for a long time. FFTs are optimized DFTs for special cases (e.g., lengths that are powers of 2).

MoonKnight has already provided a radix-4 Decimation In Frequency Cooley-Tukey scheme in Fortran. I'm below providing a radix-2 Decimation In Frequency Cooley-Tukey scheme in Matlab.

The code is an iterative one and considers the scheme in the following figure:

A recursive approach is also possible.

As you will see, the implementation calculates also the number of performed multiplications and additions and compares it with the theoretical calculations reported in How many FLOPS for FFT?.

The code is obviously much slower than the highly optimized FFTW exploited by Matlab.

Note also that the twiddle factors omegaa^((2^(p - 1) * n)) can be calculated off-line and then restored from a lookup table, but this point is skipped in the code below.

For a Matlab implementation of an iterative radix-2 Decimation In Time Cooley-Tukey scheme, please see 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); 

Your code is correct to obtain the DFT.

The function you are testing is (sin ((double) i / points * frequency * 2) which corresponds to a synoid of amplitude 1, frequency 1 and sampling frequency Fs = number of points taken.

Operating with the obtained data we have:

As you can see, the DFT coefficients are symmetric with respect to the position coefficient N / 2, so only the first N / 2 provide information. The amplitude obtained by means of the module of the real and imaginary part must be divided by N and multiplied by 2 to reconstruct it. The frequencies of the coefficients will be multiples of Fs / N by the coefficient number.

If we introduce two sinusoids, one of frequency 2 and amplitude 1.3 and another of frequency 3 and amplitude 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)); } 

The obtained data are:

Good luck.