FFT之蝶形运算

前言

写这篇博客时,是今天考完电信传输理论后的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=ej2π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=2L1B=2^{L-1}个点;

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

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

  4. 同一旋转因子对应着间隔为D=2LD=2^L点的2ML2^{M-L}个蝶形;(此点包含了第3点,但多一个信息)。

也就是说:

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

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

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

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

    2. 每个蝶形运算的两个输入数据相隔B=2L1B=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,你还这么浪,赶快滚去复习。

赞赏