C语言FFT

FFT

有具体实验要求可私聊定制

算法分析理论及仿真

实验主要通过两部分实现

$FFT$

快速傅里叶变换 (fast Fourier transform),即利用计算机计算离散傅里叶变换（$DFT)$的高效、快速计算方法的统称，简称$FFT$。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，$FFT$算法计算量的节省就越显著。

任何一个N为2的整数幂(即$N=2M$)的$DFT$，都可以通过M次分解，最后成为2点的$DFT$来计算。M次分解构成了从x(n)到X(k)的M级迭代计算，每级由N/2个蝶形组成。 例如四点$DFT$的蝶形变化图如下：

本次的实验采用了变址$FFT$运算，数据按自然顺序输入存储，然后通过“变址”运算将自然顺序转换成码位倒置顺序存储，此处采用了雷德算法实现了变址运算。

可视化

可视化部分主要依托于$EasyX$的图形库来实现，可以设置画板纬度，背景网格，线的粗细，颜色等，并依托于graphics.h库设计了绘图类，方便后续绘图调用。

实验环境

$VS2019$

实验过程与分析

模块分布

工程主要分为三个部分：主函数部分，$FFT$部分，绘图部分。FFT.c部分储存$FFT$子函数，通过直接调用子函数对导入数据进行快速傅里叶变化，graph2d.c以及graph2d.h是基于graphics.h图形库自定义的一个方便调用的二维绘图库，实现了定义并绘制坐标轴的功能，使得$FFT$波形可以更清晰地显示。主函数部分负责数据的选择和读取，本次工程采用了256点的$FFT$变化，通过$MATLAB$生成原数据文件并保存到$TXT$文件中，通过对$TXT$文件的读取，将数据读取到数组中，并进行下一步的分析处理。

代码设计与分析

C语言中并不能进行复数运算，因此需要首先定义复数运算规则，可以通过结构体实现

static struct compx plural(struct compx a, struct compx b) {
    struct compx c;
    c.real = a.real * b.real - a.imag * b.imag;
    c.imag = a.real * b.imag + a.imag * b.real;
    return(c);
}

a, b为进行运算的两个复数，c为返回的运算结果，实现了复数乘法的定义和计算

蝶形运算可以通过同址运算和变址运算实现，本次实验通过变址运算实现$FFT$算法，变址运算通过雷德算法实现

for (i = 0; i < (FFT_N-1); i++) {
        if (i < j) {    // 如果i<j,即进行变址
            t = dat[j];
            dat[j] = dat[i];
            dat[i] = t;
        }
        k = nv2;
        // 求j的下一个倒位序
        while (k <= j) // 如果k<=j,表示j的最高位为1
        {
            j = j - k; // 把最高位变成0
            k = k / 2; // k/2，比较次高位，依次类推，逐个比较，直到某个位为0
        }
        j = j + k; // 把0改为1
    }

雷德算法：

倒位序从二进制的角度来看，就是把顺序的二进制数翻转过来，以8点$DFT$为例

十进制		二进制
顺序	逆序	顺序	逆序
0	0	000	000
1	4	001	100
2	2	010	010
3	6	011	110
4	1	100	001
5	5	101	101
6	3	110	011
7	7	111	111

蝶形变化

for (l = 1; (f = f / 2) != 1; l++); 	// 蝶形级数
        for (m = 1; m <= l; m++) 		// 控制蝶形结级数
        {
            le = 2 << (m - 1);         	//le蝶形结距离，即第m级蝶形的蝶形结相距le点
            lei = le / 2; 				// 同一蝶形结中参加运算的两点的距离
            u.real = 1.0; 				// u为蝶形结运算系数，初始值为1
            u.imag = 0.0;
            w.real = cos(PI / lei); 	// w为系数商，即当前系数与前一个系数的商
            w.imag = -sin(PI / lei);
            for (j = 0; j <= lei - 1; j++) // 控制计算不同种蝶形结，即计算系数不同的蝶形结
            {
                for (i = j; i <= FFT_N - 1; i = i + le) // 控制同一蝶形结运算，即计算系数相同蝶形结
                {
                    ip = i + lei; 		// i，ip分别表示参加蝶形运算的两个节点
                    t = plural(dat[ip], u); // 蝶形运算，详见公式
                    dat[ip].real = dat[i].real - t.real;
                    dat[ip].imag = dat[i].imag - t.imag;
                    dat[i].real = dat[i].real + t.real;
                    dat[i].imag = dat[i].imag + t.imag;
                }
                u = plural(u, w); 		 // 改变系数，进行下一个蝶形运算
            }
        }

蝶形运算的第L级中：

1. 每个蝶形的输入数据相隔$B=2^{L-1}$个点；
2. 每级有B个不同的旋转因子；
3. 相同旋转因子对应的每个蝶形的输出数据相隔$D=2^L$个点；
4. 同一旋转因子对应着间隔为$D=2^L$点的$2^{M-L}$个蝶形；（此点包含了第3点，但多一个信息）。

也就是

1. 第一层循环控制级数，共$M=logN$级，从$L=1,...,M$；

2. 第二层循环控制每层的旋转因子，每层有$B=2^{L-1}$个不同的旋转因子；

3. 第三层循环控制每个旋转因子对应的蝶形运算。

   • 每个旋转因子会被使用$2^{M-L}$次；

   • 每个蝶形运算的两个输入数据相隔$B=2^{L-1}$个点，输出也是相隔B个点（蝶形：交叉平行）；

   • 而同一旋转因子对应的两个相邻的蝶形运算相隔$D=2^L$个点；

绘图头文件定义

typedef struct mypoint {
	double x;
	double y;
}Point;
class graph2d {
private:
	double height;											// 画板高度
	double width;											// 画板宽度
	Point pointlb;											// 坐标轴左下角的点
	Point pointrt;											// 坐标轴右上角的点
	int x_len;												// x轴字的宽度
	int y_len;												// y轴字的宽度
    
public:
    // 初始化
	graph2d();												// 初始化为默认值
	graph2d(double _width, double _height, Point _pointlb, Point _pointrt);		// 初始化画板宽高及网格
	~graph2d();												// 析构函数
	void waitKey(int _delay = 0);							// 等待关闭
    
private:
    // 坐标转换、绘制
	wchar_t* ctow(const char* str);							// char* to wchar_t*
	void setlen(int _len = 3);								// 设置坐标文本长度
	std::string dtos(double _num, char _axis);				// string to double
	void setGrid(Point _pointlb, Point _pointrt);			// 设置网格
	void drawGrid();										// 绘制网格
	void drawAxisX();										// 绘制X轴
	void drawAxisY();										// 绘制Y轴
	bool isBorder(Point _point);							// 是否在边界内
	int numConversion(double _num, char _axis);				// 将y轴颠倒
	Point fucCSDataToAbsCSData(Point _point);				// 方程的点转换到画板的点
	Point absCSDataToFucCSData(Point _point);				// 画板的点转换到方程的点
	void showError(std::string _err);						// 显示错误
	void drawRectangle(Point _pointlb, Point _pointrt, COLORREF _colorl, COLORREF _colorf, int _style = BS_SOLID);	//画正方形
	void setBackgroundColor(COLORREF _color = 0xEAEAEA);	// 设置画板背景颜色
	void setAxisColor();									// 设置坐标系背景颜色
	void initAxis();										// 初始化坐标轴内的信息
    
public:
    // 绘制坐标方程函数
	void plot(Point _point, COLORREF _color = RED, int _size = 3, int _type = BS_SOLID);		// 绘制点
	void plot(std::vector<Point> _point, COLORREF _color = BLACK, int _thickness = 3, int _type = PS_SOLID);		// 绘制一连串的线
	void title(std::string _str);							// 标题
	void xlabel(std::string _str);							// x轴文本注释
	void ylabel(std::string _str);							// y轴文本注释
};

画板初始化，定义长宽并绘制背景网格

void graph2d::plot(std::vector<Point> _point, COLORREF _color, int _thickness, int _type)
{
	setlinecolor(_color);
	setlinestyle(_type, _thickness);
	for (int i = 1; i < _point.size(); i++) {
		if (isBorder(_point[i - 1]) && isBorder(_point[i])) {
			line(numConversion(fucCSDataToAbsCSData(_point[i - 1]).x, 'x'), numConversion(fucCSDataToAbsCSData(_point[i - 1]).y, 'y'), numConversion(fucCSDataToAbsCSData(_point[i]).x, 'x'), numConversion(fucCSDataToAbsCSData(_point[i]).y, 'y'));
		}
		else {
			showError("line");
		}
	}
}

坐标转换

return { ((_point.x - pointlb.x) / (pointrt.x - pointlb.x) * 0.78 * width + 0.12 * width),((_point.y - pointlb.y) / (pointrt.y - pointlb.y) * 0.78 * height + 0.1 * height) };

因为在绘图界面中加入了网格和标尺作为北京，因此需要对坐标进行一定的转换来达到定位准确的目的。

析构函数

graph2d::~graph2d() {
	closegraph();
}

释放内存并等待按键按下后关闭绘图窗口。

设置标题

void graph2d::title(std::string _str) {
	RECT r = { numConversion(0,'x'), numConversion(0.88 * height,'y'), numConversion(width,'x'), numConversion(height,'y') };
	settextcolor(BLACK);
	settextstyle(int(20 * height / 590), 0, _T("宋体"));
	drawtext(ctow(_str.c_str()), &r, DT_CENTER | DT_VCENTER | DT_SINGLELINE);
}

x, y轴标题设置同理

读取数据，将$txt$文件与执行文件放在同一目录下，输入完整文件名可读取数据

cout << "请输入要读取的数据文件名\n";
char txt_name[30] = "suibainxuan";
cin >> txt_name;
fopen_s(&fp, txt_name, "r");
for (int m = 0; m < FFT_N; m++) {
    start:
    rewind(stdin);
    fscanf_s(fp, "%lf %f", &fre[m], &s[m].real);
    s[m].imag = 0;
}

转换为可视化数据

FFT(s);  //进行快速傅里叶变换
            for (i = 0; i < FFT_N; i++) //求变换后结果的模值，存入复数的实部部分
                if (i == 0) {
                    result[0] = sqrt(s[i].real * s[i].real + s[i].imag * s[i].imag) / FFT_N * 2;
                    result_new[0] = result[0];
                    cout << result[0]<<"\n";
                }
                else {
                    result[i] = sqrt(s[i].real * s[i].real + s[i].imag * s[i].imag) / FFT_N * 2;
                    result_new[i] = result[i];
                    cout << result[i]<<"\n";
                }

power[0] = 10 * log(result[0] * result[0]);
power[i] = 10 * log(result[i] * result[i]);

确认绘图边界

sort(result_new, result_new + FFT_N);
edge_downside = result_new[0];			// 上边界
edge_upside = result_new[FFT_N-1];		// 下边界
double floor(edge_downside);			// 向下取整
double ceil(edge_upside);				// 向上取整

测试参数设计

原始信号为频率为$15Hz$的正弦波与$40Hz$正弦波的叠加，一组数据不添加噪声，另一组数据添加高斯白噪声，使信噪比为$0dB$，采用$MATLAB$生成原始数据并保存到$TXT$文件内，程序读取文件内数据并进行处理，将处理所得结果分别通过数据形式和可视化图像形式输出，将数据部分保存并导入$MATLAB$中与单独利用$MATLAB$处理所得数据进行比较和平均插值运算，判断所编写程序的准确度是否满足理论需求。

参数调试

调试过程中快速傅里叶变化部分较为顺利，所得数据的精准度和分辨率均达到了预想的要求

调试重点在于绘图类，在调试绘图类的时候，需要对输入坐标进行转换，使得坐标可以正确显示在对应的网格坐标点上。

可扩展部分

在本次实验的$FFT$部分采用了结构体进行对复数运算的定义和实现，但存在一定的不便捷行，进一步可以尝试使用类定义中的运算符重载来重新定义$FFT$部分的复数运算法则，例如：

class plural{
public:
    float imag, real;					 // 虚部，实部
    plural operator*(plural const& b) {	 // 复数乘法运算符重载
		CompNum sum;
		sum.real = this->real * b.real - this->imag * b.imag;
		sum.imag = this->imag * b.imag + this->real * b.imag;
		return sum;
	}
}

但是在调试过程中出现一些问题，暂未解决。

实验结果显示与验证

$FFT$

如上图所示为分别利用自己所写程序与$MATLAB$进行快速傅里叶变化后的结果对比

再将两种方法所得结果在$MATLAB$中进行对比error = mean(1-data_C./data_mat')，得出最终结果为$4.9507e-08$，两者结果近似相同。

功率谱密度

如上图所示为分别利用自己所写程序与$MATLAB$进行快速傅里叶变化后的结果对比

再将两种方法所得结果在$MATLAB$中进行对比error = mean(1-data_C./data_mat')，结果近似为0，达到预期效果

带噪信号时频分析

• 上图分别为用C和matlab绘图生成图像