Welcome to CDJ's Page

自律更自由

FFT之蝶形运算

发布于 2021-01-20 | 标签: DSP | 19分钟 | 4033字数

前言

写这篇博客时，是今天考完电信传输理论后的5个小时，没带计算器，没找别人借，一个人抗下了所有。

用C语言实现FFT之时域抽取法的想法是在考数字信号处理前萌发的，在复习时看到书上有编程思想，然后自己又好久都没有敲代码了，手有点痒，然后就按照书上的编程思想进行coding。但是之间有许多错误，哎！用C语言这么久了，还是没有养成编写代码时注意数组越界的习惯。

昨天下午考了数字信号处理，很意外，没有一点点防备，考点不对头，有些考纲上没写，比如考纲上网络结构只写了直接型、级联型以及并联型，没有提到频率采样、线性相位，但是它！它！它！考了，而且考的好偏（也许是自己太菜了）。另外，我在很久之前就在写数字信号处理的文章/博客，感觉对考试一点帮助也没有，我考前对窗函数设计FIR滤波器很熟悉，感觉很稳，结构它考了频率采样法，我又一次死了。

考完了数字信号处理的晚上，闲的蛋痛，继续考前未完成的FFT算法，然后经过慢慢思考，发现书上也有错误，关于最内层循环的迭代方式不对，经过理解，然后修改，得出了正确结果（和MATLAB对比，当然我是用的在线octave进行计算的，我电脑打开MTALAB慢惨）。

今天，此时此刻以及前一个小时，手又有点痒了，于时打开FFT代码，继续修改，优化（其实也没怎么优化），完成了从文件中读取数据，将FFT的结果写入到文件中。

遂作感想如上，记录于此

下面继续写blog，总结。

相关文件

fft.c
fft_rawdata.txt
fft_cookeddata.txt

从fft_rawdata.txt中读取要进行fft的数据，数据格式为：

  1. 第一行一个数字，表示数据点数/长度；
  2. 此后每行有两个数字，表示实部和虚部，用Table分割；

如下：

fft的结果会保存到fft_cookeddata.txt中去，如下：

在代码中我还计算了每个fft结果的幅度和相位，但是没有保存到文件中去，而是打印到标准输出，结果如下（结果中我是以%.1f输出的，所以只显示一位小数）：

实操

准备工作

对于复数，直接用一个结构体即可，两个成员，分别表示实部和虚部，后面又定义了一个结构体，其实是一样的，只不过成员名称不一样了，用幅度和角度来表示一个复数。

typedef struct {
	float real;
	float imag;
}complex_t;

// 另外一种复数的表达，幅度加角度/相位
typedef struct {
	float amp;
	float pha;
}complex_ap_t;

时域抽取法/蝶形运算的编程思想参考了数字信号处理书上的第四章。

旋转因子的表示：

刚开始不知道WNpW_N^p怎么表示，后面一想，它不就是一个复数嘛，用欧拉公式展开就可以得到实部和虚部。

WNp=e−j2πNp=cos⁡(2πp/N)−jsin⁡(2πp/N)W_N^p = e^{-j\frac{2\pi}{N}p} = \cos(2{\pi}p/N)-j\sin(2{\pi}p/N)

因为旋转因子用的比较多，比较频繁，所以没有写成函数，而是用了一个宏定义，如下：

#define pi 	(3.1415926535)
// 旋转因子，{}语句块，其结果为最后一条语句的值
#define	WNP(N,P)	({			\
	complex_t temp;				\
	temp.real = cos(2*pi*P/N);	\
	temp.imag = -sin(2*pi*P/N);	\
	temp; })

然后还需要实现复数的基本四则运算，加、减、乘，除（代码中没有实现），还有求一个复数的幅度和相位。

倒序算法

另外，蝶形运算的初始输入的数据需要倒叙，所以需要写一个函数进行输入数据的倒叙，倒序的算法为：对原始数据的序号的二进制进行反转，得到的结果就是其倒叙后的位置（假设总共有N个数据，那么只看log2Nlog_2^N个二进制位）。

比如总共有8个数据，那么就值用看3个二进制位，第3个原始数据，经过倒叙后是第6个数据（从0开始，有第0个数据）。

3 = 011

6 = 110

而C语言不能实现数据的二进制位倒相，所以又要写一个函数（此函数来源于网络）。

//采用移位的方法使一个数的二进制位翻转后返回
unsigned int reverse_bit(unsigned int num, const unsigned int N)
{
	unsigned int ret = 0;
	int bit = 0;
	int i = 0;
	// log(8.0)/log(2.0) = 3.000,but (unsigned int)x = 2,why?
	unsigned int M = (unsigned int)(log(N+0.00001)/log(2.0));
	for (i = 0; i < M; i++)
	{
		ret <<= 1;
		bit = num & 1;
		ret = bit | ret;
		num = num / 2;
	}
	//1000 0000 0000 0000 0000 0000 0000 0000
	return ret;
}

数据的倒序：

/**
 * 进行FFT前对数据进行倒叙以方便时域抽取法
 * @Author   CofCai
 * @DateTime 2020-12-27T20:48:41+0800
 * @param    c_arr                    数据指针
 * @param    N                        数据长度
 * @return                            成功：1、失败：0
 */
int reverse_order(complex_t* c_arr, const unsigned int N)
{
	complex_t c_temp;
	/** calloc: initial by zeros  */
	char* bitmap = calloc(N, sizeof(char));
	unsigned int i, j;
	for (i = 0; i < N; ++i) {
		if (bitmap[i] == 0) {
			// 先保存第i个数据
			c_temp.real = c_arr[i].real;
			c_temp.imag = c_arr[i].imag;
			// 计算i对应的倒叙，和这个倒叙位置交换数据
			j = reverse_bit(i, N);
			c_arr[i].real = c_arr[j].real;
			c_arr[i].imag = c_arr[j].imag;
			c_arr[j].real = c_temp.real;
			c_arr[j].imag = c_temp.imag;
			// 防止后面j又和i进行倒叙交换
			*(bitmap+i) = 1;
			*(bitmap+j) = 1;
		} else {
			continue;
		}
	}
	return 1;
}

编程思想

蝶形运算要点：
在第L级中：

  1. 每个蝶形的输入数据相隔B=2L−1B=2^{L-1}个点；

  2. 每级有B个不同的旋转因子；

  3. 相同旋转因子对应的每个蝶形的输出数据相隔D=2LD=2^L个点；

  4. 同一旋转因子对应着间隔为D=2LD=2^L点的2M−L2^{M-L}个蝶形；（此点包含了第3点，但多一个信息）。

也就是说：

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

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

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

    1. 每个旋转因子会被使用2M−L2^{M-L}次；

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

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

代码实现

/**
 * author:		CofCai
 * datatime:	2020-12-27 21:00:24
 * file description:
 *	 根据数据信号处理书上第4四章FFT提到的时域抽取法编程思想，
 *	 用C语言实现DIT-FFT算法。
 *	 编写过程中遇到问题了：
 *	 	旋转因子该怎么表示?
 *	 	突然又解决了，用旋转因子不就是复数蛮，一样可以用复数表示
 *	 	W_N^k = e^{-j2pik/N} = cos(2pik/N) - j*sin(2pik/N)
 *
 * 	 似乎结果有问题
 * 	 考完了数字信号处理（老师给的考点一点也不靠谱，好多不会，awsl）
 * 	 继续敲代码，做事情前先弄懂原理，刚开始跟着书上的抄，结果不对，
 * 	 后面自己认真看书，理解蝶形运算的过程，然后发现第三层循环，也就
 * 	 是最内循环的迭代方式不对，不应该是++，而应该按照D=2^L的方式迭代
 * 	 
 * 	 蝶形运算要点：
 *	 	在第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. 第三层循环控制每个旋转因子对应的蝶形运算
 * 				1. 每个旋转因子会被使用2^{M-L}次；
 * 				2. 每个蝶形运算的两个输入数据相隔B=2^{L-1}个点，输出也是相隔B个点（蝶形：交叉平行）；
 * 				3. 而同一旋转因子对应的两个相邻的蝶形运算相隔D=2^L个点；
 *
 *
 * 	该文件的使用：
 * 		fft_rawdata.txt：此文件用于输入原始数据，
 * 						 第一行是数据个数N，此后每一行表示一个数据，分为实部和虚部，
 * 						 中间以table分割，以换行结束，如：
 * 						 3
 * 						 1.2\t2.3\n
 * 						 2.3\t78.2\n
 * 						 7.2\t12.3\n
 * 		fft_cookeddata.txt:此文件用于保存结果
 * 	
 * 	另外需要注意：
 * 		传入fft的是数据指针，因此是在原数据上直接保存结果，所以如果还需要用到原数据，那么请copy一份。
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

typedef struct {
	float real;
	float imag;
}complex_t;

// 另外一种复数的表达，幅度加角度/相位
typedef struct {
	float amp;
	float pha;
}complex_ap_t;


#define DBG_INFO 0
#if DBG_INFO
  	#define DBG_PRINTF(fmt, args...)  \
	do\
	{\
	    printf(fmt, ##args);\
	}while(0)
#else
  #define DBG_PRINTF(fmt, args...)   
#endif

#define pi 	(3.1415926535)
// 旋转因子，{}语句块，其结果为最后一条语句的值
#define	WNP(N,P)	({			\
	complex_t temp;				\
	temp.real = cos(2*pi*P/N);	\
	temp.imag = -sin(2*pi*P/N);	\
	temp; })

#define R2C(R)		({			\
	complex_t c;				\
	c.real = R;					\
	c.imag = 0;					\
	c; })

float cplx_amp(complex_t* c)
{
	return sqrt(pow(c->real, 2) + pow(c->imag, 2));
}

float cplx_angle(complex_t* c)
{
	return atan2(c->imag, c->real);
}

complex_ap_t cplx_amp_angle(complex_t* c)
{
	complex_ap_t c_ap;
	c_ap.amp = cplx_amp(c);
	c_ap.pha = cplx_angle(c);

	return c_ap;
}

// 为了减小开销，还是传入指针吧
complex_t cplx_add(complex_t* c1, complex_t* c2)
{
	complex_t temp;
	temp.real = c1->real + c2->real;
	temp.imag = c1->imag + c2->imag;
	return temp;
}

complex_t cplx_sub(complex_t* c1, complex_t* c2)
{
	complex_t temp;
	temp.real = c1->real - c2->real;
	temp.imag = c1->imag - c2->imag;
	return temp;
}

complex_t cplx_prod(complex_t* c1, complex_t* c2)
{
	complex_t temp;
	// (a1+b1j)(a2+b2j) = (a1a2-b1b2) + (a1b2 + b1a2)
	temp.real = (c1->real * c2->real) - (c1->imag * c2->imag);
	temp.imag = (c1->real * c2->imag) + (c1->imag * c2->real);
	return temp;
}

//采用移位的方法使一个数的二进制位翻转后返回
unsigned int reverse_bit(unsigned int num, const unsigned int N)
{
	unsigned int ret = 0;
	int bit = 0;
	int i = 0;
	// log(8.0)/log(2.0) = 3.000,but (unsigned int)x = 2,why?
	unsigned int M = (unsigned int)(log(N+0.00001)/log(2.0));
	for (i = 0; i < M; i++)
	{
		ret <<= 1;
		bit = num & 1;
		ret = bit | ret;
		num = num / 2;
	}
	return ret;
}

/**
 * 进行FFT前对数据进行倒叙以方便时域抽取法
 * @Author   CofCai
 * @DateTime 2020-12-27T20:48:41+0800
 * @param    c_arr                    数据指针
 * @param    N                        数据长度
 * @return                            成功：1、失败：0
 */
int reverse_order(complex_t* c_arr, const unsigned int N)
{
	complex_t c_temp;
	/** calloc: initial by zeros  */
	char* bitmap = calloc(N, sizeof(char));
	unsigned int i, j;
	for (i = 0; i < N; ++i) {
		if (bitmap[i] == 0) {
			// 先保存第i个数据
			c_temp.real = c_arr[i].real;
			c_temp.imag = c_arr[i].imag;
			// 计算i对应的倒叙，和这个倒叙位置交换数据
			j = reverse_bit(i, N);
			c_arr[i].real = c_arr[j].real;
			c_arr[i].imag = c_arr[j].imag;
			c_arr[j].real = c_temp.real;
			c_arr[j].imag = c_temp.imag;
			// 防止后面j又和i进行倒叙交换
			*(bitmap+i) = 1;
			*(bitmap+j) = 1;
		} else {
			continue;
		}
	}
	return 1;
}

/**
 * FFT算法
 * 在第L级中：
 * 			1. 每个蝶形的输入数据相隔B=2^{L-1}个点；
 * 			2. 每级有B个不同的旋转因子；
 * 			3. 相同旋转因子对应的每个蝶形的输出数据相隔D=2^L个点；
 * 		 	4. 同一旋转因子对应着间隔为D=2^L点的2^{M-L}个蝶形；（此点包含了第3点，但多一个信息）。
 * @Author   CofCai
 * @DateTime 2020-12-28T21:20:27+0800
 * @param    c_arr                    输入序列，为指针
 * @param    N                        序列长度
 * @return                            无
 */
int fft(complex_t* c_arr, const unsigned int N)
{
	//log_a^b = log_c^b / log_c^a：换底公式
	unsigned int M = (unsigned int)(log((float)N+0.00001)/log(2.0));
	DBG_PRINTF("FFT related parameter:");
	DBG_PRINTF("N = %u\tM = %d\n", N, M);
	// L:层数，最大为M
	unsigned int L, k;
	// 蝶形运算中输入的两个数据相距B个点，J用于索引不同的旋转因子
	unsigned int J, B, p;
	// 相同旋转因子的输出结果之间的间隔为：D=2^L
	unsigned int D;
	complex_t c_temp, T, WNp;
	// 倒叙
	reverse_order(c_arr, N);
	// 蝶形运算
	// 第L层
	for (L = 1; L <= M; ++L) {
		B = (unsigned int)pow(2, L-1);
		DBG_PRINTF("第 L = %d层\tB = %d个不同的旋转因子\n", L, B);
		// 第L层有B = 2^(L-1)个旋转因子，每个为W_{2^L}^J，J=0,1,...,B-1。
		// 经过变换，可得旋转因子为：W_N^p,p=Jx2^{M-L}。
		for (J = 0; J <= B-1; ++J) {
			// 第L层第J个旋转因子的指数
			p = (unsigned int)(pow(2, M-L)*J);
			DBG_PRINTF("J = %d对应的旋转因子的指数为P = %d\n", J, p);
			// 前面说过每层共有2^(L-1)个旋转因子，而每个旋转因子在此层会被利用2^{M-L}次
			WNp = WNP(N, p);
			DBG_PRINTF("旋转因子WNp: (%u, %u)-->(%1.f + %1.fj)\n", N, p, WNp.real, WNp.imag);
			// D是同一旋转因子对应的蝶形的输出结果之间的距离，而B是输入碟形数据之间的距离
			D = (unsigned int)pow(2, L);
			// 下面的循环就是计算出WNp对应的蝶形运算的结果
			for (k = J; k <= N-1; k+=D) {
				// A_{L}(J) = A_{L-1}(J) + A_{L-1}(J+B)W_N^p;
				// A_{L}(J+B) = A_{L-1}(J) - A_{L-1}(J+B)W_N^p;
				
				c_temp = cplx_prod(c_arr+k+B, &WNp);
				// DBG_PRINTF("prod:(%1.f+%1.fj) * (%1.f+%1.fj)=(%1.f+%1.fj)\n", 
				// 	(c_arr+k+B)->real, (c_arr+k+B)->imag, WNp.real, WNp.imag, c_temp.real, c_temp.imag);
				
				T = cplx_add(c_arr+k, &c_temp);
				DBG_PRINTF("AL(%u) = A(%u) + A(%u)WNp\n", k, k, k+B);
				// DBG_PRINTF("add:(%1.f+%1.fj) + (%1.f+%1.fj)=(%1.f+%1.fj)\n",
				// 	(c_arr+k)->real, (c_arr+k)->imag, c_temp.real, c_temp.imag, T.real, T.imag);
				
				*(c_arr+k+B) = cplx_sub(c_arr+k, &c_temp);
				DBG_PRINTF("AL(%u) = A(%u) - A(%u)WNp\n", k+B, k, k+B);
				// DBG_PRINTF("sub:(%1.f+%1.fj) - (%1.f+%1.fj)=(%1.f+%1.fj)\n",
				// 	(c_arr+k)->real, (c_arr+k)->imag, c_temp.real, c_temp.imag, (c_arr+k+B)->real, (c_arr+k+B)->imag);
				
				*(c_arr+k) = T;
			}
		}
		putchar('\n');
	}
	DBG_PRINTF("FFT done!!!\n\n");
	return 1;
}

void cplx_amp_angle_batch(complex_t* c, unsigned int N, complex_ap_t* c_ap)
{
	unsigned int i = 0;
	for (i = 0; i < N; ++i) {
		*(c_ap+i) = cplx_amp_angle(c+i);
	}
}

void fft_test()
{
	unsigned int i = 0, j;
	unsigned int N = 8;
	float angle = 0, abs = 0;
	complex_t c_arr[8] = {
		{1, 2},
		{3, 0},
		{5, 0},
		{7, 0},
		{1, 0},
		{3, 0},
		{5, 0},
		{7, 0}
	};
	complex_t* c_cpy = (complex_t*)calloc(N, sizeof(complex_t));
	memcpy(c_cpy, c_arr, N*sizeof(complex_t));

	complex_ap_t* c_ap = (complex_ap_t*)calloc(N, sizeof(complex_ap_t));

	fft(c_cpy, N);
	cplx_amp_angle_batch(c_cpy, N, c_ap);

	printf("i\torigin data\tfft data\tamp\tpha\n");
	printf("-\t-----------\t--------\t---\t---\n");
	for (i = 0; i < N; ++i) {
		printf("%d\t%.2f + %.2fj\t%.2f + %.2fj\t%.2f\t%.2f\n",
			i, c_arr[i].real, c_arr[i].imag, 
			(c_cpy+i)->real, (c_cpy+i)->imag,
			(c_ap+i)->amp, (c_ap+i)->pha );
	}
}


void fft_test_read_data_from_file()
{
	int i, N;
	complex_t c;
	char filename[] = "fft_rawdata.txt";
	FILE* fp = NULL;

	fp = fopen(filename, "r+");
	if (fp == NULL) {
		printf("open %s failure!\n", filename);
		exit(0);
	}
	fscanf(fp, "%d", &N);
	complex_t* c_arr = (complex_t*)calloc(N, sizeof(complex_t));
	complex_ap_t* c_ap_arr = (complex_ap_t*)calloc(N, sizeof(complex_ap_t));

	i = 0;
	while(i < N) {
		fscanf(fp, "%f\t%f", &(c_arr+i)->real, &(c_arr+i)->imag);
		printf("%.1f + %.1fj\n", c_arr[i].real, c_arr[i].imag);
		i++;
	}
	fclose(fp);

	printf("FFT begin!\n");
	fft(c_arr, N);
	printf("FFT done!\n");
	strcpy(filename, "fft_cookeddata.txt");
	fp = fopen(filename, "w+");
	if (fp == NULL) {
		printf("open %s failure!\n", filename);
		fclose(fp);
		exit(0);
	}
	i = 0;
	while(i < N) {
		c_ap_arr[i].amp = cplx_amp(c_arr+i);
		c_ap_arr[i].pha = cplx_angle(c_arr+i);
		printf("%.1f + %.1fj ---> amp=%.1f, pha=%.1f\n", 
			c_arr[i].real, c_arr[i].imag, c_ap_arr[i].amp, c_ap_arr[i].pha);
		fprintf(fp, "%f\t%f\n", c_arr[i].real, c_arr[i].imag);
		i++;
	}
	fclose(fp);
}

int main(int argc, char const *argv[])
{
	fft_test();
	// fft_test_read_data_from_file();

	return 0;
}

fft_test()函数的执行结果：

总结

fft函数传入的是数据指针，因此会修改原数据，如果还需要使用到原始数据，那么请copy一份。

总结个锤子，后天考移动通信，cdj，你还这么浪，赶快滚去复习。

赞赏
# DSP
Linux高性能服务器设计 上一篇
CMake的使用及pcre2 下一篇
已运行:
浏览数: 次 | 访客数: 人
0%
召唤看板娘