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).
among them,$W_N = e^{j\frac{2\pi}{N}}$.$X(k)$ is$x(n)$Discrete Fourier transform of$x$$($$n$$)$ .
The matrix equation can be used to see the transformation process of DFT more clearly:
$X =/begin{pmatrix} X(0)/X(1)/x(2)\\vdots/X(N1)\\end{pmatrix}$;$W =/begin{pmatrix} 1 & 1 & 1 &/cdots & 1/1 & W_N^1 & W_N^2 &/cdots & W_N^{N1}/1 & W_N^2 & W_N^4 &/cdots & W_N^{2(N1)}\\vdots &/vdots &/vdots &/ddots &/vdots/1 & W_N^{N1} & W_N^{2(N1 )} &/cdots & W_N^{(N1)(N1)}\\end{pmatrix}$;$x =/begin{pmatrix} x(0)/x(1)/x(2)\\vdots/x(N1)\\end{pmatrix}$
It can be seen that the length is$N$Finite sequence of$N$$x(n)$ , its discrete Fourier transform$X(k)$ still a length of$N$A finite sequence of$N.$It 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
1.2.1 Nature foreshadowing
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:
It can be seen that$x(n)$The DFT of$x$$($$n$$)$ becomes the combination of the DFT of the evennumbered term and the DFT of the oddnumbered 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$x$$($$n$$)$ 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 evennumbered DFT, use$H(k)$ replaces the oddnumbered DFT, then the finishing formulas (2) and (3) are:
among them
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 (arbr, aibi);}
Cpx operator * (Cpx a, Cpx b) { return Cpx (ar * brai * 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 fixedpoint 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:
Functions used when reading parameters:
//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 (arbr, aibi);}
Cpx operator * (Cpx a, Cpx b) { return Cpx (ar * brai * 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.