 ## Fast Fourier Transform (FFT) and C Implementation Using the Octave GNU Tool

### 2. Introduction

##### These transformations are quite time-consuming and to understand this, we will describe (1) in more detail: ##### There are algorithms that simplify the calculation of DFT (IDFT), such algorithms are called Fast Fourier Transform (FFT). The idea is that the existing input array is split into shorter subsequences, for which the Fourier Transform is calculated. For example, we divide the original sequence of length N into two subsequences with length N/2 each. Using (4) and (6) we obtain for each subsequence: ##### In the future array length will be used: ### 3. Decimation In Time FFT. Algorithm from Cooley-Tukey

##### Then (1) can be rewritten using (10) and (11): ##### Often this basic operation, described by formula (15), is represented as a graphical butterfly (flow graph), shown in Figure 3-1. Figure 3-1: Butterfly of the Decimation In Time FFT

##### We divide the already obtained subsequences into new ones consisting of even and odd elements. Continuing this process, we will come to a two-points FFT. See the Figure 3-2. Figure 3-2: Division the Input Sequence x(n) for the Decimation In Time FFT

##### Now using (15) and starting with the two-points Fourier Transform, we can build a diagram Decimation In Time FFT. See the Figure 3-3. Figure 3-3: 16 points of the Decimation In Time FFT

1. ##### Bit Reversal Permutation. The Figure 3-2 and Figure 3-3 are shown that the input array must be reordered: x(0), x(8), x(4), x(12),…x(7), x(15). This corresponds to the bit reversal of the sequence indices. Below in Table 3-1 is an example for a 4-bit index for a 16-point FFT. Table 3-1: Bit Reversal Example for the 4 Bit Number (16 Point FFT)

### 4. Decimation In Frequency FFT. Algorithm from Sande-Tukey

##### Figure 4-1: Butterfly of the Decimation In Frequency FFT

##### Using the same formulas (27) and (28) to subsequences we will get four subsequences with length N/4 each. Continuing this process, we will come to subsequences with length 1. These obtained values are the spectrum of the original sequence x(n), mixed up in accordance with the bit reversal. See the Figure 4-2. Figure 4-2: Division the Input Sequence x(n) for the Decimation In Frequency FFT. The spectrum of the obtained subsequences is shown on the figure

##### Now we can build a diagram of the Decimation In Frequency FFT. See the Figure 4-3. Figure 4-3: 16 points of the Decimation In Frequency FFT

### 5. Realization of the Inverse Fourier Transform Using the FFT Algorithm

##### Let’s introduce notation for complex-conjugate operation with *: ##### Then using (29), we can rewrite IDFT (2) as: ##### If x(n) is real signal then (30) can be written in more simple form: ### 7. FIR Filter. Fast Convolution using FFT

##### The traditional method of calculating the FIR filter is performed using the formula (35). ##### Fast Convolution использует FFT and the N value is selected as a power of two: ##### complex multiplications. Multiplication of the spectra is also used, this is N more complex multiplications. Each complex multiplication corresponds to four multiplications of real numbers. Then: ##### Figure 7-1: Number of multiplications for fast convolution (red) and FIR filter (blue) using (35)

### 8. C-Source Code fft_float.c/h

1. ##### Defines and constants:

#define FFT_LOG2_POINTS_NMB 5u – FFT size definition. Hier is 2^5=32 FFT points

flt_complex_nmb_t flt_complex_exp_table[FFT_LOG2_POINTS_NMB] – FFT rotation exponent factors. The table is prepared from the Octave GNU Tool, see the FloatExp2file( )

2. ##### Private Function

/*
Example: 001 => 100
Input
uint16 previous_index – previous bit reversal value
uint16 number – FFT number = 2^M
Return
uint16 – next bit reversal value
*/

3. ##### Export Functions

/*
Complex Multiplication
(a1+jb1)(a2+jb2)= a1 x a2 – b1 x b2 + j(a1 x b2 + a2 x b1)
Input
flt_complex_nmb_t* ptr_complex_nmb1 – pointer on the operand_1
flt_complex_nmb_t* ptr_complex_nmb2 – pointer on the operand_2
Return
flt_complex_nmb_t – (operand_1 x operand_2)
*/
flt_complex_nmb_t flt_complex_mult(flt_complex_nmb_t* complex_nmb1, flt_complex_nmb_t* complex_nmb2)

/*
a1+jb1 + a2+jb2 = a1+a2 + j(b1+b2)
Input
flt_complex_nmb_t* ptr_complex_nmb1 – pointer on the operand_1
flt_complex_nmb_t* ptr_complex_nmb2 – pointer on the operand_2
Return
flt_complex_nmb_t – (operand_1 + operand_2)
*/

/*
Complex Subtraction
(a1+jb1) – (a2+jb2) = a1-a2 + j(b1-b2)
Input
flt_complex_nmb_t* ptr_complex_nmb1 – pointer on the operand_1
flt_complex_nmb_t* ptr_complex_nmb2 – pointer on the operand_2
Return
flt_complex_nmb_t – (operand_1 – operand_2)
*/
flt_complex_nmb_t flt_complex_sub(flt_complex_nmb_t* complex_nmb1, flt_complex_nmb_t* complex_nmb2)

/*
Complex Conjugate
conjugate(a + jb) = a – jb
Input
flt_complex_nmb_t* ptr_complex_nmb1 – pointer on the operand
Return
flt_complex_nmb_t – conjugate(operand)
*/
flt_complex_nmb_t flt_complex_conjugate(flt_complex_nmb_t* complex_nmb)

/*
Fast Fourier Transform (FFT), Decimation In Time, Cooley-Tukey
Input
flt_complex_nmb_t* ptr_input_array – input array
flt_complex_nmb_t* ptr_exp_array – pointer on the complex exponential function
uint16 log2_size – FFT size: N = 2^log2_size
Return
void
FFT result in the ptr_input_array
*/
void flt_fft_decimation_in_time(flt_complex_nmb_t* ptr_input_array, flt_complex_nmb_t* ptr_exp_array, uint16 log2_size)

/*
Fast Fourier Transform (FFT), Decimation In Frequency, Sande–Tukey
Input
flt_complex_nmb_t* ptr_input_array – input array
flt_complex_nmb_t* ptr_exp_array – pointer on the complex exponential function
uint16 log2_size – FFT size: N = 2^log2_size
Return
void
FFT result in the ptr_input_array
*/
void flt_fft_decimation_in_frequency(flt_complex_nmb_t* ptr_input_array, flt_complex_nmb_t* ptr_exp_array, uint16 log2_size)

/*
FIR Filter using FFT: Fast Convolution
Input
float* fir_coeff_ptr – pointer on the filter coefficients
uint16 fir_order – FIR order or coefficients number
sint16* input_ptr – pointer on the input signal array
uint16 input_length – input array length
sint16 output_ptr – pointer on the output array. Note output array shall be cleared: output = {0}
uint16 log2_size – FFT size: N = 2^log2_size
Return
uint8 – TRUE is OK: FIR filter is done / FALSE is NOT OK: FFT size is not suitable for the fir order FIR output in the output_ptr array
*/
uint8 fir_filter_using_fft(float* fir_coeff_ptr, uint16 fir_order, sint16* input_ptr, uint16 input_length, sint16* output_ptr)

int main(void) – only test

### 9. C-Source Code fft_integer.c/h

1. ##### Defines and constants

#define FFT_LOG2_POINTS_NMB 5u – FFT size: 2^5=32 FFT points
#define EXP_LOG2_FACTOR 15u – Factor: 2^15=32768 => 32767*exp(-j*2*pi*n/N)
int_complex_nmb_t int_complex_exp_table[FFT_POINTS_NMB_DIV2] – FFT: N/2 exponents

2. ##### Private Functions:

/*
Example: 001 => 100
Input
uint16 previous_index – previous bit reversal value
uint16 number – FFT number = 2^FFT_LOG2_POINTS_NMB
Return
uint16 – next bit reversal value
*/

/*
Number Limit: sint32 to sint16
Input
sint32 number – input value
Return
sint16: SINT16_MIN_NEGATIVE_NMB <= output_number <= SINT16_MAX_POSITIVE_NMB
*/
sint16 number_limit(sint32 number)

3. ##### Export Functions

/*
Complex Multiplication
(a1+jb1)(a2+jb2)= a1 x a2 – b1 x b2 + j(a1 x b2 + a2 x b1)
Multiplication shall be done with the A*exp(-j*2*pi*n/N)
and is used the round and 2^EXP_LOG2_FACTOR normalization
Input
int_complex_nmb_t* ptr_complex_nmb1 – pointer on the operand_1
int_complex_nmb_t* ptr_complex_nmb2 – pointer on the operand_2
Return
int_complex_nmb_t – (operand_1 x operand_2)
*/
int_complex_nmb_t int_complex_mult(int_complex_nmb_t* complex_nmb1, int_complex_nmb_t* complex_nmb2)

/*
a1+jb1 + a2+jb2 = a1+a2 + j(b1+b2)
Input
int_complex_nmb_t* ptr_complex_nmb1 – pointer on the operand_1
int_complex_nmb_t* ptr_complex_nmb2 – pointer on the operand_2
Return
int_complex_nmb_t – (operand_1 + operand_2)
*/

/*
Complex Subtraction
(a1+jb1) – (a2+jb2) = a1-a2 + j(b1-b2)
Input
int_complex_nmb_t* ptr_complex_nmb1 – pointer on the operand_1
int_complex_nmb_t* ptr_complex_nmb2 – pointer on the operand_2
Return
int_complex_nmb_t – (operand_1 – operand_2)
*/
int_complex_nmb_t int_complex_sub(int_complex_nmb_t* complex_nmb1, int_complex_nmb_t* complex_nmb2)

/*
Complex Conjugate
conjugate(a + jb) = a – jb
The function is used for the Inverse FFT
Input
int_complex_nmb_t* ptr_complex_nmb1 – pointer on the operand
Return
int_complex_nmb_t – conjugate(operand)
*/
int_complex_nmb_t int_complex_conjugate(int_complex_nmb_t* complex_nmb)

/*
Fast Fourier Transform (FFT), Decimation In Time, Cooley-Tukey
Input
int_complex_nmb_t* ptr_input_array – input array
int_complex_nmb_t* ptr_exp_array – pointer on the complex exponential function
uint16 log2_size – FFT size: N = 2^log2_size
Return
void
FFT result in the ptr_input_array
*/
void int_fft_decimation_in_time(int_complex_nmb_t* ptr_input_array, int_complex_nmb_t* ptr_exp_array, uint16 log2_size)

/*
Fast Fourier Transform (FFT), Decimation In Frequency, Sande–Tukey
Input
int_complex_nmb_t* ptr_input_array – input array
int_complex_nmb_t* ptr_exp_array – pointer on the complex exponential function
uint16 log2_size – FFT size: N = 2^log2_size
Return
void
FFT result in the ptr_input_array
*/
void int_fft_decimation_in_frequency(int_complex_nmb_t* ptr_input_array, int_complex_nmb_t* ptr_exp_array, uint16 log2_size)

int main(void) – only test

### 10. Octave GNU Tool Script FFT_Support.m

1. ##### Constants:

fft_power = 5; % 2^5=32 FFT Points

fir_coeff_original = […]; % Fast Convolution support
F_sample_in_Hz = 12e3; % Fast Convolution support

2. ##### Functions:

% FFT Float Complex Exponent Rotation Table to File
% Support C-Code
% See the fft_float.c
% Array flt_complex_exp_table
% Input parameters:
% fftPower – FFT power: N=2^fftPower
% FileNameString – output file name
function FloatExp2file(fftPower, FileNameString)

% FFT Integer Complex Exponent Table to File
% Support C-Code
% See the fft_integer.c
% Array int_complex_exp_table
% Input parameters:
% fftPower – FFT power: N=2^fftPower
% FileNameString – output file name
function IntegerExp2file(fftPower, IntegerFactor, FileNameString)

% FFT Float Test Signal Table to File
% Support C-Code
% See the fft_float.c
% Array flt_test_signal1/flt_test_signal2
% Input parameters:
% fftPower – FFT power: N=2^fftPower
% FileNameString – output file name
% Output:
% txt file with the signal data
% signal_sum = 10sin(2*pi*n/N) + 10sin(2*pi*n/N)
function signal_sum = FloatTestSignal2file(fftPower, FileNameString)

% FFT Integer Test Signal Table to File
% Support C-Code
% See the fft_integer.c
% Array int_test_signal1/int_test_signal2
% Input parameters:
% fftPower – FFT power: N=2^fftPower
% FileNameString – output file name
% Output:
% txt file with the signal data
% signal_sum = 512*sin(2*pi*n/N) + 512*sin(2*pi*n/N)
function signal_sum = IntegerTestSignal2file(fftPower, FileNameString)

%————————–
% Fast Convolution Support
% FIR Low Pass Filter
%————————–

% Integer Test Signal to File
% Support C-Code
% Input parameters:
% SignalAmpl – amplitude of the test signal
% SignalFreqInHz – frequency of the test signal in Hz
% SignalLength – table length
% FsampleInHz – sampling frequency in Hz
% FileNameString – output file name
% Output parameters:
% TestSignal – test signal vector
% txt file
function TestSignal = IntegerTestSignals2file(SignalAmpl, SignalFreqInHz, SignalLength, FsampleInHz, FileNameString)

% Integer Table to File
% Support C-Code
% Input parameters:
% ParameterVector – data array
% FileNameString – output file name
function IntegerParamVector2file(ParameterVector, FileNameString)

% Float Table to File
% Support C-Code
% Input parameters:
% ParameterVector – data array
% FileNameString – output file name
function FloatParamVector2file(ParameterVector, FileNameString)