FFT之蝶形运算

CofCai

已于 2022-03-20 11:01:54 修改

阅读量1.1w

点赞数 12

分类专栏：信号处理文章标签： fft c语言

于 2020-12-30 16:05:27 首次发布

本文链接：https://blog.csdn.net/weixin_44360592/article/details/111992566

版权

信号处理专栏收录该内容

4 篇文章 2 订阅

订阅专栏

文章目录

Author🚹:CofCai

Email✉️:cai.dongjun@nexuslink.cn

QQ😙:1664866311

CSDN Page🏠

前言

写这篇博客时，是今天考完电信传输理论后的5个小时，没带计算器。

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

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

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

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

遂作感想如上，记录于此

下面继续写blog，总结。

实操

准备工作

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

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

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

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

旋转因子的表示：

刚开始不知道 $W_N^p$ 怎么表示，后面一想，它不就是一个复数嘛，用欧拉公式展开就可以得到实部和虚部。
$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个数据，那么只看 $log_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;
	}
	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级中：

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

也就是说：

第一层循环控制级数，共 $M = l o g N$ 级，从 $L = 1, . . ., M$ ；
第二层循环控制每层的旋转因子，每层有 $B=2^{L-1}$ 个不同的旋转因子；
第三层循环控制每个旋转因子对应的蝶形运算。
1. 每个旋转因子会被使用 $2^{M-L}$ 次；
2. 每个蝶形运算的两个输入数据相隔 $B=2^{L-1}$ 个点，输出也是相隔B个点（蝶形：交叉平行）；
3. 而同一旋转因子对应的两个相邻的蝶形运算相隔 $D=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()函数的执行结果：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ldTCAMos-1609315289573)(imgassets/fft_result_with_amp_pha.jpg)]

总结

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

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

fft_test()函数的执行结果：

在这里插入图片描述

CofCai

关注

12
点赞
踩
82

收藏

觉得还不错? 一键收藏
5
评论
FFT之蝶形运算

文章目录前言相关文件实操准备工作倒序算法编程思想代码实现总结Author????:CofCaiEmail✉️:cai.dongjun@nexuslink.cnQQ????:1664866311personal Page????CSDN Page????前言写这篇博客时，是今天考完电信传输理论后的5个小时，没带计算器，没找别人借，一个人抗下了所有。用C语言实现FFT之时域抽取法的想法是在考数字信号处理前萌发的，在复习时看到书上有编程思想，然后自己又好久都没有敲代码了，手有点痒，然后就按照书
复制链接

扫一扫