One: FFT principle
1.1 DFT calculations
It is suggested that you can follow the formula when you have time to write it, which is helpful to understand ~. The discrete Fourier series (DFS) transform over a period is defined as the discrete Fourier Transform (DFT).
The WN = e – j2 PI NW_N = e ^ {- j \ frac {2 \ PI} {N}} WN = e – jN2 PI. X(k)X(k)X(k) is the discrete Fourier transform of X(n) X(n) X(n).
The transformation process of DFT can be seen more clearly with matrix equation:
X = (X (0) (1) (2) ⋮ X X X (N – 1)) X = \ the begin (0) \ \ {pmatrix} X X (1) (2) \ \ \ \ X \ vdots \ \ X (N – 1) \ \ X = {pmatrix} \ end ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎛ X (0) (1) (2) ⋮ X X X (N – 1) ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎞; W = (111. 11 wn1wn2. WNN – 11 wn2wn4.. WN2 (N – 1) ⋮ ⋮ ⋮ ⋱ ⋮ WNN 1-1 WN2 (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)} \\ W = {pmatrix} \ end ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎛ 111 ⋮ 11 wn1wn2 ⋮ WNN – 11 wn2wn4 ⋮ WN2 (N – 1)….. ⋱. 1 WNN – 1 WN2 (N – 1) ⋮ WN (N – 1) (N – 1) ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎞; X = (x (0) (1) (2) ⋮ x x x (N – 1)) x = \ the begin (0) \ \ {pmatrix} x x (1) (2) \ \ \ \ x \ vdots \ \ x (N – 1) \ \ X = {pmatrix} \ end ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎛ x (0) (1) (2) ⋮ x x x (N – 1) ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎞
It can be seen that the finite length sequence x(n)x(n)x(n) x(n) of length NNN, its discrete Fourier transform x(k) x(k) x(k) is still a finite length sequence of length NNN. It can be seen from (1) that the time complexity is O(N2)O(N^2)O(N2). If N=1024N =1024N =1024 points, 1048576(more than one million) complex complex multiplications are needed. The DFT is so computative that it has an optimized version: the Fast Fourier Transform (FFT).
1.2 FFT computation
1.2.1 Property bedding
Because the coefficient WNnk=e−j2πNnkW_N^{nk} = E ^{-j\frac{2\ PI}{N}nk}WNnk=e−jN2πnk is a periodic function, we can use its properties to improve the algorithm and improve the computational efficiency.
-
Properties: WNk + N2 = – WNkW_N ^ {k + \ frac {N} {2}} = – W_N ^ kWNk + 2 N = – WNk (symmetry)
-
WNnk=W1nkNW_N^{nk} = W_1^{\frac{nk}{N}}WNnk=W1Nnk
Here, the above two properties are mainly used to decompose the DFT operation of large points with length of N points into several DFT of small points successively. Because the DFT is proportional to N2N^2N2, NNN is also small.
1.2.2 Basis 2FFT extracted by time (N points)
Assuming that the number N of FFT is 2 to the integral power (base 2), the sequence is first divided into two groups, one group is the even term and the other is the odd term, and then the following transformation is carried out, and the derivation is as follows:
It can be seen that finding the DFT of x(n)x(n) becomes finding the DFT of its even-numbered and odd-numbered DFT, but notice that this only computs the first half of the DFT, and the second half is obtained by the following property:
By property one:
So now we have the full DFT of x(n)x(n)x(n), and that’s the whole idea of FFT, but let’s use a butterfly diagram to make it a little bit more intuitive.
1.3 Butterfly signal flow diagram
G(k)G(k)G(k) is used to replace even term DFT, and H(k)H(k)H(k) H(k)H(k) is used to replace odd term DFT, then the collation equations (2) and (3) are:
Among them
It can be seen from (4) and (5) that we can divide a string of time domain data into even and odd parts to calculate G(K)G(K)G(K) and H(K) H(K) H(K) H(K). Similarly, we can divide even parts into even parts and odd parts to calculate until there are only two data left, and then recursively calculate FFT results. See the following classic N-point butterfly diagram for the specific straight point flow:
Two: FFT C++ implementation
#include <iostream> // FFT algorithm implementation, base 2 time extraction
#include <vector>
#include <ctime>
using namespace std;
const double PI = acos(- 1); / / the value of PI
struct Cpx// Define a complex structure and a complex operator {
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(a.r + b.r, a.i + b.i); }
Cpx operator - (Cpx a, Cpx b) { return Cpx(a.r - b.r, a.i - b.i); }
Cpx operator * (Cpx a, Cpx b) { return Cpx(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r); }
void fft(vector<Cpx>& a, int lim, int opt)
{
if (lim == 1) return;
vector<Cpx> a0(lim >> 1).a1(lim >> 1); // Initialize half of the size to store even and odd numbers
for (int i = 0; i < lim; i += 2)
a0[i >> 1] = a[i], a1[i >> 1] = a[i + 1]; // Divide into even and odd parts
fft(a0, lim >> 1, opt); // Calculate the even part recursively
fft(a1, lim >> 1, opt); // Calculate the even part recursively
Cpx wn(cos(2 * PI / lim), opt * -sin(2 * PI / lim)); / / equal to WN
Cpx w(1.0);
for (int k = 0; k < (lim >> 1); k++) // See graph 1
{
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 the butterfly figure 2, small 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(a)
{
int opt = 1; // 1 is FFT, -1 is IFFT
vector<Cpx> a(16); // This is fixed at 16 points and can be changed
for (int i = 0; i < 16; i++) // Randomly generate 16 numbers as data to be processed
{
Cpx c = Cpx(cos(0.2 * PI * i), 0);
a[i] = c;
}
if (1 == opt)
fft(a, 16, opt); // The a array becomes the value after FFT
else if (- 1 == opt)
{
fft(a, 16, opt); // The a array becomes the value after IFFT
for (int i = 0; i < 512; i++) a[i].r /= 512, a[i].i /= - 512.;// IFFT is divided by the length
}
else;
return 0;
}
Copy the codeThree: MATLAB and C++ mixed programming
In engineering sometimes in order to make data processing faster or support some fixed-point operations, and choose to use C/C++ to deal with some processing steps, in fact, the general engineering with MATLAB processing speed is enough, mixed programming is also a review of C++.
MATLAB and C++ mixed programming is divided into MATLAB call C++ and C++ call MATLAB, here we discuss the former. The hybrid programming of MATLAB and C++ 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 compiler configuration of MATLAB to the interruption point debugging of the final program jump to VS, I encountered many difficulties in the whole process of mixed programming. There are many materials available on the Internet but it is also messy. Here is a summary of the steps I have done from the beginning to the end.
I am using MATLAB2019b and VS2019, before using MATLAB2016, and then download what 2019 support files, modify the registry, etc., did not work for a long time, simply change MATLAB2019b directly.
(2) run mex-setup C++ and mbuild-setup C++ in MATLAB. If not, the current version of MATLAB does not support the current version of Visual Studio. It is not recommended to reduce the version of VS, there will be compatibility issues.
③ no need to create a project, directly create a xx. CPP file according to the mex interface definition to write a C++ program (specific procedures discussed later). I’ve been messing around with configuration issues in VS for a long time, like linking extern libraries and so on, but I don’t feel like I need to create a project in the end, so I don’t need to configure these external link libraries? Can you just write xx.cpp? (I’m not sure, it might work)
④ After the program is written in MATLAB run mex-g xx. CPP, if xx. CPP program written in accordance with the specification, it will be mex success, generate XX.mexw64 and xx.pdb files; If MEX fails, modify the code according to the warning returned by MATLAB. Note that in order to be able to debug breakpoints later in VS2019, -g is required.
⑤ Write the corresponding test program in the MATLAB script, set breakpoints to run and stop at xx() function.
⑥ Open xx. CPP file with VS2019, find ‘Add to process’ in the ‘Debug’ column, select ‘native’, and add MATLAB to process. Set breakpoints where you want to stop.
⑦ If MATLAB continues to run, it will enter the corresponding breakpoint in VS2019. (The last two steps may not go in, in fact, I can sometimes go in and sometimes not, there is no good solution for the moment)
3.2 Interface Usage
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++ program.
The most important header files and interface main functions are as follows, written in a fixed way.
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
Copy the codeNRHS: Number of input parameters
PRHS []: Array of Pointers to input parameters
NLHS: indicates the number of output parameters
PLHS [] : array of Pointers to output parameters
Note that input and output are transmitted in the form of Pointers, which can be understood as MATLAB puts its parameters at a certain address, and then C++ reads the corresponding length of data at the corresponding address according to the length of this parameter, thus completing the transmission process of parameters. Instead, pass it back at the end. Here are some commonly used MEX functions:
Functions used to read arguments:
// Complex single-value reading
double Nr1 = *mxGetPr(prhs[0]); // Read the real part of the first argument
double Ni2 = *mxGetPr(prhs[1]); // Read the imaginary part of the second argument
Copy the code// Address read
double* Pr1 = mxGetPr(prhs[0]); // Read the real part of the first argument
double* Pi2 = mxGetPi(prhs[0]); // Read the virtual address of the first argument
Copy the code// Matrix dimension read
int M = mxGetM(prhs[2]); // Read the number of rows for the third argument
int N = mxGetN(prhs[2]); // Read the number of columns for the third argument
Copy the codeTo be added
Functions used to print parameters:
// Output the 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 codeTo be added
3.3 FFT MATLAB/C++ hybrid implementation
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 structure and a complex operator {
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(a.r + b.r, a.i + b.i); }
Cpx operator - (Cpx a, Cpx b) { return Cpx(a.r - b.r, a.i - b.i); }
Cpx operator * (Cpx a, Cpx b) { return Cpx(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r); }
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]; // Divide into even and odd parts
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
{
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]); // Enter the number of rows in the matrix
int N = mxGetN(prhs[0]); // Enter the number of matrix columns
double* xpr = mxGetPr(prhs[0]); // Input matrix real part pointer
double* xpi = mxGetPi(prhs[0]); // Enter the imaginary part of the matrix pointer
int lim = *mxGetPr(prhs[1]); Lim = N, M = 1
int opt = *mxGetPr(prhs[2]); // Select 1 for FFT and -1 for IFFT
plhs[0] = mxCreateDoubleMatrix(M, N, mxCOMPLEX); // Output matrix creation (critical)
double* ypr = mxGetPr(plhs[0]); // Output the real part of the matrix pointer
double* ypi = mxGetPi(plhs[0]); // Outputs the virtual part of the matrix pointer
std::vector<Cpx> a(lim); // Use vector to store data
for (int i = 0; i < lim; i++) // Input vector is passed in
{
a[i].r = xpr[i];
a[i].i = xpi[i];
}
if (1 == opt)
fft(a, lim, opt); // The a array changes to the value after FFT
else if (- 1 == opt)
{
fft(a, lim, opt); // The a array changes to the value after IFFT
for (int i = 0; i < lim; i++) a[i].r /= lim, a[i].i /= lim;// IFFT is divided by the length
}
else;
for (int i = 0; i < lim; i++) // The output vector comes out
{
ypr[i] = a[i].r;
ypi[i] = a[i].i;
}
return;
}
Copy the codeThen write the following program in MATLAB script:
clear all
mex fftxx.cpp -g
a = randn(1.16) + 1i * randn(1.16); % Randomly generate 16 complex data
fftsize = 16;
b = fftxx(a, fftsize, 1) % is passed into C++ for FFT processing
b1 = fft(a,fftsize) % MATLAB system functions for FFT processing
c = fftxx(b, fftsize, - 1) % is passed into C++ for IFFT processing
c1 = ifft(b, fftsize) % MATLAB system functions for IFFT processing
Copy the codeFinally, run the. M program, in the MATLAB command line window, it can be seen that b and B1, c and c1 output results are completely consistent.
Four: reference materials
www.bilibili.com/video/BV1Y7…