（C++实现fft的方案）Matlab转C的方案总结_matlab转化fft函数为c语言

导言

 经过了很长一段时间的实验，总结下matlab转C给Qt使用的方案总结，这里的实验对象为matlab的fft函数。

方案对比

FFT翻译方案对比

C++的fft普通实现

//快速傅里叶变换fft
   QVector<Complex> fft(const QVector<double> &real, const QVector<double> &imag) {
   	/*
   		real：向量实部
   		imag：向量虚部
   		这里仅实现对向量的快速傅里叶变换：Y = fft(X)
   		计算公式可见matlab关于fft的介绍

   		Y(k) = n∑(j=1) X(j)*Wn^(j-1)(k-1)
   		X(j) = n∑(k=1)Y(k)*Wn^(j-1)(k-1)
   		其中 Wn = e^((-2πi)/n) = cos(-2π/n)+sin(-2π/n)i
   	*/
   	QVector<Complex> out(real.size(), {0.0, 0.0});
   	double fixed_factor = (-2 * PI) / real.size();
   	for (int u = 0; u < real.size(); u++) {
   		double uxf = u * fixed_factor;
   		for (int x = 0; x < real.size(); x++) {
   			//X(j)*Wn^(j-1)(k-1)
   			double power = uxf * x;
   			double temp = imag[x] * sin(power);
   			out[u].real += real[x] * cos(power) - temp;
   			out[u].imag += real[x] * sin(power) + temp;
   		}
   	}

   	return out;
   }
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

蝶形运算的fft

#include <list>
#include <cmath>
#include <string>
#include <vector>
#include <iostream>

#define PI acos(-1)
using namespace std;

//定义复数结构
struct Complex
{
   double imag;
   double real;
};

//复数乘法函数
Complex ComplexMulti(Complex One, Complex Two)
{
   //Temp用来存储复数相乘的结果
   Complex Temp;
   Temp.imag = One.imag * Two.real + One.real * Two.imag;
   Temp.real = One.real * Two.real - One.imag * Two.imag;
   return Temp;
}

//此处的length为原序列的长度     将输入的任意数组填充，以满足快速傅里叶变换
void Data_Expand(double *input, vector<double> &output, const int length)
{
   int i = 1, flag = 1;
   while (i < length)
   {
   	i = i * 2;
   }

   //获取符合快速傅里叶变换的长度
   int length_Best = i;

   if (length_Best != length)
   {
   	flag = 0;
   }

   if (flag)
   {
   	//把原序列填充到Vector中
   	for (int j = 0; j < length; ++j)
   	{
   		output.push_back(input[j]);
   	}
   }

   else
   {
   	//把原序列填充到Vector中
   	for (int j = 0; j < length; ++j)
   	{
   		output.push_back(input[j]);
   	}
   	//把需要拓展的部分进行填充
   	for (int j = 0; j < length_Best - length; j++)
   	{
   		output.push_back(0);
   	}
   }
}

//此处的length为填充后序列的长度//将输入的实数数组转换为复数数组
void Real_Complex(vector<double> &input, Complex *output, const int length)
{
   for (int i = 0; i < length; i++)
   {
   	output[i].real = input[i];
   	output[i].imag = 0;
   }
}

//FFT变换函数//快速傅里叶变换函数
//input: input array
//StoreResult: Complex array of output
//length: the length of input array
void FFT(Complex *input, Complex *StoreResult, const int length)
{

   //获取序列长度在二进制下的位数
   int BitNum = log2(length);

   //存放每个索引的二进制数，重复使用，每次用完需清零
   list<int> Bit;

   //遍历序列的每个索引
   for (int i = 0; i < length; ++i)
   {
   	//flag用来将索引转化为二进制
   	//index用来存放倒叙索引
   	//flag2用来将二值化的索引倒序
   	int flag = 1, index = 0, flag2 = 0;

   	//遍历某个索引二进制下的每一位
   	for (int j = 0; j < BitNum; ++j)
   	{
   		//十进制转化为长度一致的二进制数，&位运算符作用于位，并逐位执行操作
   		//从最低位（右侧）开始比较
   		//example:
   		//a = 011, b1 = 001, b2 = 010 ,b3 = 100
   		//a & b1 = 001, a & b2 = 010, a & b3 = 000
   		int x = (i & flag) > 0 ? 1 : 0;

   		//把x从Bit的前端压进
   		Bit.push_front(x);

   		//等价于flag = flag << 1，把flag的值左移一位的值赋给flag
   		flag <<= 1;
   	}

   	//将原数组的索引倒序，通过it访问Bit的每一位
   	for (auto it : Bit)
   	{
   		//example:其相当于把二进制数从左到右设为2^0，2^1，2^2
   		//Bit = 011
   		//1: index = 0 + 0 * pow(2, 0) = 0
   		//2: index = 0 + 1 * pow(2, 1) = 2
   		//3: index = 2 + 1 * pow(2, 2) = 6 = 110
   		index += it * pow(2, flag2++);
   	}

   	//把Data[index]从DataTemp的后端压进，其实现了原序列的数据的位置调换
   	//变换前：f(0), f(1), f(2), f(3), f(4), f(5), f(6), f(7)
   	//变换后：f(0), f(4), f(2), f(6), f(1), f(5), f(3), f(7)
   	StoreResult[i].real = input[index].real;
   	StoreResult[i].imag = input[index].imag;

   	//清空Bit
   	Bit.clear();
   }

   //需要运算的级数
   int Level = log2(length);
   Complex Temp, up, down;

   //进行蝶形运算
   for (int i = 1; i <= Level; i++)
   {
   	//定义旋转因子
   	Complex Factor;

   	//没有交叉的蝶形结的距离(不要去想索引)
   	//其距离表示的是两个彼此不交叉的蝶型结在数组内的距离
   	int BowknotDis = 2 << (i - 1);

   	//同一蝶形计算中两个数字的距离（旋转因子的个数）
   	//此处距离也表示的是两个数据在数组内的距离(不要去想索引)
   	int CalDis = BowknotDis / 2;

   	for (int j = 0; j < CalDis; j++)
   	{
   		//每一级蝶形运算中有CalDis个不同旋转因子
   		//计算每一级（i）所需要的旋转因子
   		Factor.real = cos(2 * PI / pow(2, i) * j);
   		Factor.imag = -sin(2 * PI / pow(2, i) * j);

   		//遍历每一级的每个结
   		for (int k = j; k < length - 1; k += BowknotDis)
   		{
   			//StoreResult[k]表示蝶形运算左上的元素
   			//StoreResult[k + CalDis]表示蝶形运算左下的元素
   			//Temp表示蝶形运算右侧输出结构的后半部分
   			Temp = ComplexMulti(Factor, StoreResult[k + CalDis]);

   			//up表示蝶形运算右上的元素
   			up.imag = StoreResult[k].imag + Temp.imag;
   			up.real = StoreResult[k].real + Temp.real;

   			//down表示蝶形运算右下的元素
   			down.imag = StoreResult[k].imag - Temp.imag;
   			down.real = StoreResult[k].real - Temp.real;

   			//将蝶形运算输出的结果装载入StoreResult
   			StoreResult[k] = up;
   			StoreResult[k + CalDis] = down;
   		}
   	}
   }
}

int main()
{
   //获取填充后的Data;
   vector<double> Data;
   //进行FFT的数据，此处默认长度为2的幂次方
   //double Datapre[] = {1, 2, 6, 4, 48, 6, 7, 8, 58, 10, 11, 96, 13, 2, 75, 16};
   double Datapre[] = { 100, 2, 56, 4, 48, 6, 45, 8, 58, 10 };
   //数据扩充
   Data_Expand(Datapre, Data, 10);

   //将输入数组转化为复数数组
   Complex array1D[16];
   //存储傅里叶变换的结果
   Complex Result[16];
   //变为复数形式
   Real_Complex(Data, array1D, 16);

   //快速傅里叶变换
   FFT(array1D, Result, 16);

   //基于范围的循环，利用auto自动判断数组内元素的数据类型
   for (auto data : Result)
   {
   	//输出傅里叶变换后的幅值
   	cout << data.real << "\t" << data.imag << endl;
   }

   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
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
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214

注意

该文章仅个人学习使用，欢迎大家一起交流学习