# FFT principle and detailed introduction of C++ and MATLAB mixed programming

## One: FFT principle

### 1.1 DFT calculation

It is recommended that you have time to follow the formula to write it down, which is helpful for understanding~ . The discrete Fourier series (DFS) transform in a period is defined as the discrete Fourier transform (DFT).

$\begin{cases} X(k) =/sum_{n=0}^{N-1}x(n)W_N^{kn}, & 0/le k/le {N-1}/x(n ) =/frac{1}{N}/sum_{k=0}^{N-1}X(k)W_N^{-kn}, & 0/le n/le {N-1}\\end {cases}$

among them,$W_N = e^{-j\frac{2\pi}{N}}$.$X(k)$ is$x(n)$Discrete Fourier transform of .

The matrix equation can be used to see the transformation process of DFT more clearly:

\begin{aligned} X = W/cdot x\tag{1}/end{aligned}

$X =/begin{pmatrix} X(0)/X(1)/x(2)\\vdots/X(N-1)\\end{pmatrix}$;$W =/begin{pmatrix} 1 & 1 & 1 &/cdots & 1/1 & W_N^1 & W_N^2 &/cdots & W_N^{N-1}/1 & W_N^2 & W_N^4 &/cdots & W_N^{2(N-1)}\\vdots &/vdots &/vdots &/ddots &/vdots/1 & W_N^{N-1} & W_N^{2(N-1 )} &/cdots & W_N^{(N-1)(N-1)}\\end{pmatrix}$;$x =/begin{pmatrix} x(0)/x(1)/x(2)\\vdots/x(N-1)\\end{pmatrix}$

It can be seen that the length is$N$Finite sequence of$x(n)$ , its discrete Fourier transform$X(k)$ still a length of$N$A finite sequence ofIt can be seen from (1) that the time complexity is$O(N^2)$if$N = 1024$ points, 1,048,576 (more than one million) complex number multiplications are required. The amount of calculation of DFT is really too large, so there is a later optimized version: Fast Fourier Transform (FFT).

### 1.2 FFT calculation

Due to the coefficient$W_N^{nk} = e^{-j\frac{2\pi}{N}nk}$ is a periodic function, and its properties can be used to improve the algorithm and increase the computational efficiency.

• Nature 1:$W_N^{k +/frac{N}{2}} = -W_N^k$ (symmetry)

• Nature 2:$W_N^{nk} = W_1^{\frac{nk}{N}}$ (Divide by a constant)

Here, the above two properties are mainly used to decompose the DFT operation of a large number of points with a length of N points into a number of DFTs of small points in turn. Because the amount of DFT calculation is proportional to$N^2$ ,$N$ small calculations are also small.

#### 1.2.2 Radix 2 FFT (N points) decimated by time

Assuming that the number of FFT points N is a power of 2 (base 2), first divide the sequence into two groups, one for even items and one for odd items, and then perform the following transformations, which are derived as follows:

\begin{aligned} X(k) &=/sum_{n=0}^{N-1}x(n)W_N^{nk}/&=/sum_{n=0 is an even number}^{N- 2}x(n)W_{N}^{nk} +/sum_{n=1 is an odd number}^{N-2}x(n)W_{N}^{nk}/&=/sum_{r =0}^{\frac{N}{2}-1}x(2r)W_{N}^{2rk} +/sum_{r=0}^{\frac{N}{2}-1}x (2r+1)W_{N}^{(2r+1)k}/&=/sum_{r=0}^{\frac{N}{2}-1}x(2r)W_{N} ^{2rk} + W_N^k/sum_{r=0}^{\frac{N}{2}-1}x(2r+1)W_{N}^{2rk} (by nature two/downarrow)//&=/sum_{r=0}^{\frac{N}{2}-1}x(2r)W_{\frac{N}{2}}^{rk} + W_N^k/sum_{r =0}^{\frac{N}{2}-1}x(2r+1)W_{\frac{N}{2}}^{rk}, 0/le k/le/frac{N}{ 2}-1/tag{2}\\end{aligned}

It can be seen that$x(n)$The DFT of becomes the combination of the DFT of the even-numbered term and the DFT of the odd-numbered term. Note that this only calculates the first half of the DFT value, and the second half is obtained by the following properties:

By nature one:



In this way, we have calculated the complete$x(n)$The DFT value of is actually the core idea of FFT. Next, we use butterfly diagrams to make the above calculation steps more intuitive.

### 1.3 Butterfly signal flow diagram

use$G(k)$ instead of even-numbered DFT, use$H(k)$ replaces the odd-numbered DFT, then the finishing formulas (2) and (3) are:

$\begin{cases} X(k) = G(k) + W_N^k H(k), & 0/le k/le/frac{N}{2}-1\X(k +/frac{N }{2}) = G(k)-W_N^k H(k), & 0/le k/le/frac{N}{2}-1/tag{4}\\end{cases}$

among them

$\begin{cases} G(k) =/sum_{r=0}^{\frac{N}{2}-1}x(2r)W_{\frac{N}{2}}^{rk}//H(k) =/sum_{r=0}^{\frac{N}{2}-1}x(2r+1)W_{\frac{N}{2}}^{rk}/tag{ 5}\\end{cases}$

As can be seen from (4) and (5), we can divide a string of time domain data into even and odd parts to calculate$G(K)$ and$H(k)$ , the even part can also be divided into even part and odd part again, until there are only two data left at the end, and then recursively calculate the FFT result. For the specific intuitive point of the process, see the classic N point below Butterfly chart:

## Two: C++ implementation of FFT

# include  <iostream>  //fft algorithm implementation, base 2 time extraction
# include  <vector>
# include  <ctime>
using  namespace std;

const  double PI = acos ( -1 ); //pi value

struct  Cpx//Define a complex number structure and complex number arithmetic
{
double r, i;
Cpx (): r ( 0 ), i ( 0 ) {)
Cpx ( double _r, double _i): r (_r), i (_i) {}
};
Cpx operator + (Cpx a, Cpx b) { return  Cpx (ar + br, ai + bi);}
Cpx operator- (Cpx a, Cpx b) { return  Cpx (ar-br, ai-bi);}
Cpx operator * (Cpx a, Cpx b) { return  Cpx (ar * br-ai * bi, ar * bi + ai * br);}

void  fft (vector<Cpx>& a, int lim, int opt)
{
if (lim == 1 ) return ;
vector<Cpx> a0 (lim >> 1 ) , a1 (lim >> 1 ) ;//initialize half Size, storing even and odd parts
for ( int i = 0 ; i <lim; i += 2 )
a0[i >> 1 ] = a[i], a1[i >> 1 ] = a[i + 1 ]; //divided into even part and odd part

fft (a0, lim >> 1 , opt); //recursively calculate the even part
fft (a1, lim >> 1 , opt); //recursively calculate the even part

Cpx wn (cos( 2 * PI/lim), opt * -sin( 2 * PI/lim)) ; //equals to WN
Cpx w ( 1 , 0 ) ;
for ( int k = 0 ; k <(lim >> 1 ); k++) //See butterfly diagram 1 operation process
{
a[k] = a0[k] + w * a1[k];
a[k + (lim >> 1 )] = a0[k]-w * a1[k];
w = w * wn;
}

//for (int k = 0; k <(lim >> 1); k++)//See butterfly diagram 2, a little optimization, one less multiplication
//{
//Cpx t = w * a1[k];
//a[k] = a0[k] + t;
//a[k + (lim >> 1)] = a0[k]-t;
//w = w * wn;
//}

}

int  main ()
{
int opt = 1 ;//1 is FFT, -1 is IFFT
vector<Cpx> a ( 16 ) ;//here is fixed at 16 points, which can be changed
for ( int i = 0 ; i < 16 ; i++)//Randomly generate 16 numbers as the data to be processed
{
Cpx c = Cpx ( cos ( 0.2 * PI * i), 0 );
a[i] = c;
}

if ( 1 == opt)
fft (a, 16 , opt); //a array becomes the value after FFT
else  if ( -1 == opt)
{
fft (a, 16 , opt); //a array becomes the value after IFFT
for ( int i = 0 ; i < 512 ; i++) a[i].r/= 512 , a[i].i/= - 512 ; //IFFT must be divided by length
}
else ;

return  0 ;
}
Copy code

## Three: MATLAB and C++ mixed programming

Sometimes in engineering, in order to make data processing faster or support certain fixed-point calculations, some processing steps are selected to be processed by C/C++. In fact, the processing speed of MATLAB is sufficient for general engineering, and mixed programming is also considered as it is. Review C++.

The mixed programming of MATLAB and C++ is divided into calling C++ in MATLAB and calling MATLAB in C++. Here we are discussing the former. MATLAB and C++ mixed programming is not simply to write the two languages together, but to follow an interface specification, which is discussed in 3.2.

### 3.1 Mixed programming steps

From the MATLAB compiler configuration to the final program jump to the breakpoint debugging in VS, I encountered a lot of difficulties in the whole process of mixed programming. There are many materials that can be found on the Internet but it is also messy. Here is a summary of what I have been from the beginning to The last steps.

I used MATLAB2019b and VS2019. I used MATLAB2016 before. Then I downloaded what support files for 2019, modified the registry, etc. After a long time, I didn't get it right, so I just changed to MATLAB2019b.

Run mex -setup C++ and mbuild -setup C++ in MATLAB. If it is unsuccessful, it means that the current version of MATLAB does not support the current version of Visual Studio. It is recommended to upgrade the MATLAB version. It is not recommended to lower the version of VS, as there will be compatibility problems.

No need to create a project, directly create a xx.cpp file and write a C++ program according to the mex interface definition (the specific program will be discussed later). I created the project for a long time to troubleshoot the configuration problems in VS, such as linking the extern library, etc., but I feel that in the end there is no need to create a project, so there is no need to configure these external link libraries? Just write the xx.cpp file directly? (I m not so sure, it may be useful)

After the program is written, run mex -g xx.cpp in MATLAB. If the xx.cpp program is written in compliance with the specification, the mex will succeed and generate the xx.mexw64 and xx.pdb files; if the mex fails, it will be returned according to MATLAB Warning to modify the code. Note that in order to be able to enter the breakpoint debugging in VS2019 later, you must add -g.

Write the corresponding test program in the MATLAB script, set the breakpoint to run and stop at the xx() function.

Open the xx.cpp file with VS2019, find "Add to Process" in the "Debug" column, go in and select "Local", and then add MATLAB to the process. Set a breakpoint where you want to stop.

When MATLAB continues to run, it will enter the corresponding breakpoint in VS2019. (The last two steps may not be able to get in, in fact, sometimes I can get in and sometimes I can t, and there is no good solution for the time being)

### 3.2 Interface use

The mex file is the bridge between the .m file in MATLAB and the .cpp file in VS. The quality of the mex interface is related to whether our MATLAB data can run correctly in C++ programs.

The most important header files and interface main functions are as follows, and the writing method is fixed.

# include  "mex.h"

void  mexFunction ( int nlhs, mxArray plhs * [], int nrhs, const mxArray * prhs [])
Copy the code

nrhs: the number of input parameters

prhs[]: pointer array of input parameters

nlhs: the number of output parameters

plhs[]: pointer array of output parameters

Note that input and output are transmitted in the form of pointers. It can be understood that MATLAB puts its parameters at a certain address, and then C++ reads the data of the corresponding length at the corresponding address according to the length of this parameter, and it is completed The process of passing parameters. Instead, pass it back in the end. A few commonly used mex functions are summarized below:

//Read the complex single value
double Nr1 = * mxGetPr (prhs[ 0 ]); //Read the real part of the first parameter
double Ni2 = * mxGetPr (prhs[ 1 ]); //Read the second parameter
Copy code of imaginary part
//Read the address
double * Pr1 = mxGetPr (prhs[ 0 ]); //Read the real address of the first parameter
double * Pi2 = mxGetPi (prhs[ 0 ]); //Read the first parameter
Copy code of imaginary address
//read matrix dimensions
int M = mxGetM (prhs[ 2 ]); //read the number of rows of the third parameter
int N = mxGetN (prhs[ 2 ]); //read the number of columns of the third parameter
Copy code

Functions used when outputting parameters:

//Output complex matrix
plhs[ 0 ] = mxCreateDoubleMatrix (M, N, mxCOMPLEX); //Create M*N complex matrix
double * outPr = mxGetPr (plhs[ 0 ]);
double * outPi = mxGetPi (plhs[ 0 ] );
copy the code

### 3.3 MATLAB/C++ hybrid implementation of FFT

First rewrite the FFT code in Chapter 2 into the following form using mex interface, and save it as fftxx.cpp.

# include  "mex.h"
# include  <vector>
# include  <ctime>

const  double PI = acos ( -1 ); //pi

struct  Cpx//Define a complex number structure and complex number arithmetic
{
double r, i;
Cpx (): r ( 0 ), i ( 0 ) {}
Cpx ( double _r, double _i): r (_r), i (_i) {}
};
Cpx operator + (Cpx a, Cpx b) { return  Cpx (ar + br, ai + bi);}
Cpx operator- (Cpx a, Cpx b) { return  Cpx (ar-br, ai-bi);}
Cpx operator * (Cpx a, Cpx b) { return  Cpx (ar * br-ai * bi, ar * bi + ai * br);}

void  fft (std::vector<Cpx>& a, int lim, int opt)
{
if (lim == 1 ) return ;
std::vector<Cpx> a0 (lim >> 1 ) , a1 (lim >> 1 ) ;
for ( int i = 0 ; i <lim; i += 2 )
a0[i >> 1 ] = a[i], a1[i >> 1 ] = a[i + 1 ]; //divided into even part and odd part

fft (a0, lim >> 1 , opt);
fft (a1, lim >> 1 , opt);

Cpx wn (cos( 2 * PI/lim), opt * -sin( 2 * PI/lim)) ;
Cpx w ( 1 , 0 ) ;
for ( int k = 0 ; k <(lim >> 1 ); k++ ) //Butterfly operation process
{
a[k] = a0[k] + w * a1[k];
a[k + (lim >> 1 )] = a0[k]-w * a1[k];
w = w * wn;
}
}

void  mexFunction ( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])  //mex main function
{
int M = mxGetM (prhs[ 0 ]);//input matrix row number
int N = mxGetN ( prhs[ 0 ]);//Enter the number of matrix columns
double * xpr = mxGetPr (prhs[ 0 ]);//Enter the pointer to the real part of the matrix
double * xpi = mxGetPi (prhs[ 0 ]);//Enter the pointer to the imaginary part of the matrix
int lim = * mxGetPr (prhs[ 1 ]);//Input parameters, length, here is a row vector, so lim = N, M = 1
int opt = * mxGetPr (prhs[ 2 ]); //Input parameters, select, 1 is FFT, -1 is IFFT

plhs[ 0 ] = mxCreateDoubleMatrix (M, N, mxCOMPLEX); //output matrix creation (important)
double * ypr = mxGetPr (plhs[ 0 ]); //pointer to the real part of the output matrix
double * ypi = mxGetPi (plhs[ 0 ]); //pointer to the imaginary part of the output matrix

std::vector<Cpx> a (lim) ; //use vector to store data
for ( int i = 0 ; i <lim; i++) //input vector
{
a[i].r = xpr[i];
a[i].i = xpi[i];
}

if ( 1 == opt)
fft (a, lim, opt); //a array becomes the value after FFT
else  if ( -1 == opt)
{
fft (a, lim, opt); //a array becomes the value after IFFT
for ( int i = 0 ; i <lim; i++) a[i].r/= lim, a[i].i/= lim; //IFFT must be divided by length
}
else ;

for ( int i = 0 ; i <lim; i++) //output vector out
{
ypr[i] = a[i].r;
ypi[i] = a[i].i;
}

return ;
}
Copy code

Then write the following program in the MATLAB script:

clear all
mex fftxx.cpp -g
a = randn ( 1 , 16 ) + 1 i * randn ( 1 , 16 ); % randomly generate 16 complex numbers
fftsize = 16 ;
b = fftxx(a, fftsize, 1 )    % Pass into C++ for FFT processing
b1 = fft(a,fftsize)          % MATLAB system function for FFT processing

c = fftxx(b, fftsize, -1 )   % Pass into C++ for IFFT processing
c1 = ifft(b, fftsize)        % MATLAB system function for IFFT processing

Copy code

Finally, run the .m program, and you can see that the output results of b and b1, c and c1 are exactly the same in the MATLAB command line window.

## Four: Reference materials

www.bilibili.com/video/BV1Y7...