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

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).

{X(k)= n=0N 1x(n)WNkn,0 k N 1x(n)=1N k=0N 1X(k)WN kn,0 n N 1\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,WN=e j2 NW_N = e^{-j\frac{2\pi}{N}}.X(k)X(k) isx(n)x(n)Discrete Fourier transform of .

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

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

X=(X(0)X(1)x(2) X(N 1))X =/begin{pmatrix} X(0)/X(1)/x(2)\\vdots/X(N-1)\\end{pmatrix};W=(111...11WN1WN2...WNN 11WN2WN4...WN2(N 1) 1WNN 1WN2(N 1)...WN(N 1)(N 1))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=(x(0)x(1)x(2) x(N 1))x =/begin{pmatrix} x(0)/x(1)/x(2)\\vdots/x(N-1)\\end{pmatrix}

It can be seen that the length isNNFinite sequence ofx(n)x(n) , its discrete Fourier transformX(k)X(k) still a length ofNNA finite sequence ofIt can be seen from (1) that the time complexity isO(N2)O(N^2)ifN=1024N = 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 coefficientWNnk=e j2 NnkW_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:WNk+N2= WNkW_N^{k +/frac{N}{2}} = -W_N^k (symmetry)

  • Nature 2:WNnk=W1nkNW_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 toN2N^2 ,NN 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:

X(k)= n=0N 1x(n)WNnk= n=0Is evenN 2x(n)WNnk+ n=1OddN 2x(n)WNnk= r=0N2 1x(2r)WN2rk+ r=0N2 1x(2r+1)WN(2r+1)k= r=0N2 1x(2r)WN2rk+WNk r=0N2 1x(2r+1)WN2rk(By nature two )= r=0N2 1x(2r)WN2rk+WNk r=0N2 1x(2r+1)WN2rk,0 k N2 1(2)\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 thatx(n)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 completex(n)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

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

{X(k)=G(k)+WNkH(k),0 k N2 1X(k+N2)=G(k) WNkH(k),0 k N2 1(4)\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

{G(k)= r=0N2 1x(2r)WN2rkH(k)= r=0N2 1x(2r+1)WN2rk(5)\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 calculateG(K)G(K) andH(k)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:

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

To be added

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

To be added

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...