FFT原理 & C++实现简单FFT代码_c++库函数 实现fft

傅里叶变换的意义

为什么我们要用正弦曲线来代替原来的曲线呢？

用正余弦来表示原信号会更加简单，因为正余弦拥有其他信号所不具备的性质：正弦曲线保真度。一个正弦曲线信号输入后，输出的仍是正弦曲线，只有幅度和相位可能发生变化，但是频率和波的形状仍是一样的，且只有正弦曲线才拥有这样的性质，正因如此我们才不用方波或三角波来表示。

傅立叶变换的物理意义在哪里？

傅立叶原理表明：任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加。而根据该原理创立的傅立叶变换算法利用直接测量到的原始信号，以累加方式来计算该信号中不同正弦波信号的频率、振幅和相位。当然这是从数学的角度去看傅立叶变换。

所以，最前面的时域信号在经过傅立叶变换的分解之后，变为了不同正弦波信号的叠加，我们再去分析这些正弦波的频率，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。

傅立叶变换提供给我们这种换一个角度看问题的工具，看问题的角度不同了，问题也许就迎刃而解！

FFT原理

计算机只能处理离散的数值信号，我们的最终目的是运用计算机来处理信号的。所以对于离散信号的变换只有离散傅立叶变换（DFT）才能被适用，对于计算机来说只有离散的和有限长度的数据才能被处理，对于其它的变换类型只有在数学演算中才能用到，在计算机面前我们只能用DFT方法，我们要讨论的FFT也只不过是DFT的一种快速的算法。

DFT的运算过程是这样的：

可见，在计算机上进行的DFT，使用的输入值是数字示波器经过ADC后采集到的采样值，也就是时域的信号值，输入采样点的数量决定了转换的计算规模。变换后的频谱输出包含同样数量的采样点，但是其中有一半的值是冗余的，通常不会显示在频谱中，所以真正有用的信息是N/2+1个点。

FFT的过程大大简化了在计算机中进行DFT的过程，简单来说，如果原来计算DFT的复杂度是NN次运算（N代表输入采样点的数量），进行FFT的运算复杂度是Nlg10(N)，因此，计算一个1,000采样点的DFT，使用FFT算法只需要计算3,000次，而常规的DFT算法需要计算1,000,000次！

典型的时域2分裂算法图示如下：

具体原理可以找一本数字信号处理的书籍来看。

频谱泄露：

所谓频谱泄露，就是信号频谱中各谱线之间相互干扰，使测量的结果偏离实际值，同时在真实谱线的两侧的其它频率点上出现一些幅值较小的假谱。产生频谱泄露的主要原因是采样频率和原始信号频率不同步，造成周期的采样信号的相位在始端和终端不连续。简单来说就是因为计算机的FFT运算能力有限，只能处理有限点数的FFT，所以在截取时域的周期信号时，没有能够截取整数倍的周期。信号分析时不可能取无限大的样本。只要有截断不同步就会有泄露。

避免频谱泄露的方法除了尽量使采集速率与信号频率同步之外，还可以采用适当的窗函数。

C++代码实现

只是参考大佬封装的代码雏形，后期优化地方很多，仅供参考

My_FFT.h

#pragma once
#include <cmath>
#include <vector>


#define pi 3.1415926
using namespace std;

struct My_Complex {
	double real;
	double imag; //a:real b:Imagine
};


class My_FFT {

public:
	My_FFT();
	~My_FFT();

	
	vector<My_Complex> FFT(vector<double>& res);
	//返回对应频率和幅值
	vector<pair<double, double>> abs_FFT();

	//复数求和
	My_Complex comp_plus(My_Complex u, My_Complex v);

	//复数相乘
	My_Complex comp_times(My_Complex u, My_Complex v);
	//复数相减
	My_Complex comp_minus(My_Complex u, My_Complex v);
	//序数重排
	void rev();
	My_Complex Wn(double A, double B);
	void fft1(int l, int r, int len);
public:
	int N;						//Num of Sample Nodes
	double fs=512;					//Sample frequency
	vector<My_Complex> x, u, W;
	double f0;
};

My_FFT::My_FFT() {

}
My_FFT::~My_FFT() {

}


vector<My_Complex> My_FFT::FFT(vector<double>& res) {

	this->N = int(res.size());
	this->fs = fs;
	this->f0 = fs / N;

	My_Complex temp;

	for (int i = 0; i < N; i++) {

		u.push_back(temp);
		temp.real = res[i];
		temp.imag = 0;
		
		x.push_back(temp);
		
		W.push_back(Wn(N, i));//e^-j((2*pi*i)/N)
	}

	this->rev(); //reverse the original order 
	this->fft1(0, N - 1, N);


	return x;
}

vector<pair<double, double>> My_FFT::abs_FFT() {

	vector<pair<double, double>> ret;

	for (int i = 0; i < N ; i++) {
		double u, f;
		u = sqrt(x[i].real * x[i].real + x[i].imag * x[i].imag)*2/N;
		f = double(i) * f0;
		ret.push_back({ f ,u });
	}

	return ret;
}




//复数求和
My_Complex My_FFT::comp_plus(My_Complex u, My_Complex v) {
	My_Complex e;
	e.real = u.real + v.real;
	e.imag = u.imag + v.imag;
	return e;
}

//复数相乘
My_Complex My_FFT::comp_times(My_Complex u, My_Complex v) {
	My_Complex e;
	e.real = u.real * v.real - u.imag * v.imag;
	e.imag = u.real * v.imag + u.imag * v.real;
	return e;
}

//复数相减
My_Complex My_FFT::comp_minus(My_Complex u, My_Complex v) {
	My_Complex e;
	e.real = u.real - v.real; e.imag = u.imag - v.imag;
	return e;
}


//序数重排
void My_FFT::rev() { //reverse the sequence in bit order.  
	int len = log2(N);
	int idd, ret, bit;
	for (int i = 0; i < N; i++) {
		idd = i; ret = 0;
		for (int j = 0; j < len; j++) {
			ret = ret << 1;
			bit = idd & 1;
			idd = idd >> 1;
			ret = bit | ret;
		}
		u[i] = x[ret];
	}
	for (int i = 0; i < N; i++)
		x[i] = u[i];
}

My_Complex My_FFT::Wn(double A, double B) {
	My_Complex u;
	u.real = cos((2 * pi / A) * B); 
	u.imag = -sin((2 * pi / A) * B);
	return u;
}

void My_FFT::fft1(int l, int r, int len) {
	My_Complex tmp;
	int n = len / 2;
	if (len < 2) return; //Level 1 

	fft1(l, l + n - 1, n); //even seg process 
	fft1(l + n, r, n);	//odd seg process

	for (int i = l; i <= r; i++) {
		if (i < l + n) {
			u[i] = comp_plus(x[i], comp_times(W[(i - l) * (N / len)], x[i + n]));
		}
		else {
			x[i] = comp_minus(x[i - n], comp_times(x[i], W[(i - n - l) * (N / len)]));
			x[i - n] = u[i - n];
		}
	}
	return;
}


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164

main.cpp

#include <iostream>
#include <fstream>
#include <sstream>
#include <cmath>
#include <string>
#include <vector>


#include "My_FFT.h"

using namespace std;




int main() {

	fstream  inFile1("./emulated_data.txt", ios::in);//创建一个fstream文件流对象
	vector<double> temp;
	vector<double> vx, vy, vz, p;
	string line;
	string field;
	while (getline(inFile1, line))//getline(inFile, line)表示按行读取CSV文件中的数据
	{
		line += " ";
		istringstream sin(line); //将整行字符串line读入到字符串流sin中
		while (getline(sin, field, '\t'))
		{
			temp.push_back(atof(field.c_str()));
		}
		vx.push_back(temp[0]);
		vy.push_back(temp[1]);
		vz.push_back(temp[2]);
		p.push_back(temp[3]);
		temp.clear();
	}

	My_FFT vxfft;
	vxfft.FFT(vx);
	vector<pair<double, double>> t;
	t=vxfft.abs_FFT();

	for (int i = 0; i < 512; i++) {
		cout << t[i].second << endl;//利用对象outfile往文件中写入数据
	}

	ofstream outfile("FD.txt", ios::trunc);//data.txt里面的内容会被清空

	for (int i = 0; i < 512; i++) {
		outfile << t[i].second << endl;//利用对象outfile往文件中写入数据
	}

	outfile.close();

	ofstream outfile1("TD.txt", ios::trunc);//data.txt里面的内容会被清空

	for (int i = 0; i < 512; i++) {


		outfile1 << vx[i] << endl;//利用对象outfile往文件中写入数据
		//the real Ampltude of the signal is the Amp of fft sequence times 2 divide N.
	}

	outfile1.close();

	return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67

参考Couriersix
参考https://blog.csdn.net/shenziheng1/article/details/52891807