【安富莱DSP教程】第30章 复数FFT的实现
特别说明：完整45期数字信号处理教程，原创高性能示波器代码全开源地址：链接第30章 复数FFT的实现
本章主要讲解复数FFT的浮点和定点Q31，Q15的实现。
本章节使用的复数FFT函数来自ARM官方库的TransformFunctions部分
30.1 复数FFT
30.2 复数FFT-基2算法
30.3 复数FFT-基4算法
30.4 总结
30.1 复数FFT
30.1.1 描述
当前复数FFT函数支持三种数据类型，分别是浮点，Q31和Q15。这些FFT函数有一个共同的特点，就是用于输入信号的缓冲，在转化结束后用来存储输出结果。这样做的好处是节省了RAM空间，不需要为输入和输出结果分别设置缓存。由于是复数FFT，所以输入和输出缓存要存储实部和虚部。存储顺序如下：{real, imag, real, imag,………………} ，在使用中切记不要搞错。
30.1.2 浮点
浮点复数FFT使用了一个混合基数算法，通过多个基8与单个基2或基4算法实现。根据需要，该算法支持的长度和每个长度使用不同的旋转因子表。
浮点复数FFT使用了标准的FFT定义，FFT正变换的输出结果会被放大fftLen倍数，计算FFT逆变换的时候会缩小到1/fftLen。这样就与教科书中的定义一致了。
定义好的旋转因子和位反转表已经在头文件arm_const_structs.h中定义好了，调用浮点FFT函数arm_cfft_f32时，包含相应的头文件即可。比如：
arm_cfft_f32(arm_cfft_sR_f32_len64, pSrc, 1, 1)
上式就是计算一个64点的FFT逆变换包括位反转。数据结构arm_cfft_sR_f32_len64可以认为是常数，计算的过程中是不能修改的。同样是这种数据结构还能用于混合基的FFT正变换和逆变换。
早期发布的浮点复数FFT函数版本包含基2和基4两种方法实现的，但是不推荐大家再使用了。现在全部用arm_cfft_f32代替了。
30.1.3 定点Q31和Q15
定点库提供了基2和基4两种算法，基2算法支持的数据长度，基4算法支持的数据长度。一般情况下，建议使用基4算法，基4算法比基2算法执行速度要快一些。
为了防止计算结果溢出，定点FFT每个蝶形运算的结果都要做放缩处理。对于基2算法，每次蝶形运算的结果要做0.5倍的放缩。对于基4算法，每次蝶形运算的结果要做0.25倍的放缩。FFT逆变换也要做相同的处理。相对于标准教科书式的FFT定义，定点FFT的计算结果放缩了1/fftLen（数据）倍。定点FFT的逆变也要做放缩处理，但是跟教科书式的FFT定义是相符的。
每个FFT变换都需要一个单独的结构体，但结构体中的旋转因子和位反转表可以被重新使用。
每个数据类型都有一个相关的初始化函数，初始化函数主要完成如下操作：
l 初始化结构体成员。
l 初始化旋转因子和位反转表指针。
使用初始化函数是可选的。尽管如此，如果使用了初始化函数，那么结构体不能放在const data段，如果要放在const data段，应当按照如下方法进行初始化：
*arm_cfft_radix2_instance_q31 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor};
*arm_cfft_radix2_instance_q15 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor};
*arm_cfft_radix4_instance_q31 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor};
*arm_cfft_radix4_instance_q15 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor};
*arm_cfft_instance_f32 S = {fftLen, pTwiddle, pBitRevTable, bitRevLength};
Q15和Q31 FFT使用了一个比较大的旋转因数和位反转表。这个表定义了一个最大长度的变换，表的子集可以用于短的变化。
30.1.4 arm_cfft_f32
函数定义如下：
void arm_cfft_f32(
const arm_cfft_instance_f32 * S,
float32_t * p1,
uint8_t ifftFlag,
uint8_t bitReverseFlag)
参数定义：
*S points to an instance of the floating-point CFFT structure.
*p1 points to the complex data buffer of size <code>2*fftLen</code>. Processing
occurs in-place.
ifftFlag flag that selects forward (ifftFlag=0) or inverse (ifftFlag=1) transform.
bitReverseFlag flag that enables (bitReverseFlag=1) or disables (bitReverseFlag=0) bit
reversal of output.
注意事项：
结构const arm_cfft_instance_f32的定义如下（在文件arm_math.h文件）：
typedef struct
{
uint16_t fftLen;
const float32_t *pTwiddle;
const uint16_t *pBitRevTable;
uint16_t bitRevLength;
} arm_cfft_instance_f32;
下面通过在开发板上运行这个函数并计算幅频相应，然后再与Matlab计算的结果做对比。
#define TEST_LENGTH_SAMPLES 2048
/* 输入和输出缓冲 */
static float32_t testOutput;
static float32_t testInput_f32_10khz;
/* 变量 */
uint32_t fftSize = 1024;
uint32_t ifftFlag = 0;
uint32_t doBitReverse = 1;
/*
*********************************************************************************************************
* 函 数 名: arm_cfft_f32_app
* 功能说明: 调用函数arm_cfft_f32_app计算幅频。
* 形 参：无
* 返 回 值: 无
*********************************************************************************************************
*/
void arm_cfft_f32_app(void)
{
uint16_t i;
/* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
for(i=0; i<1024; i++)
{
/* 虚部全部置零 */
testInput_f32_10khz = 0;
/* 50Hz正弦波，采样率1KHz ,作为实部 */
testInput_f32_10khz = arm_sin_f32(2*3.1415926f*50*i/1000);
}
/* CFFT变换 */
arm_cfft_f32(&arm_cfft_sR_f32_len1024, testInput_f32_10khz, ifftFlag, doBitReverse);
/* 求解模值*/
arm_cmplx_mag_f32(testInput_f32_10khz, testOutput, fftSize);
/* 串口打印求解的模值 */
for(i=0; i<1024; i++)
{
printf("%f\r\n", testOutput);
}
}
运行如上函数可以通过串口打印出计算的模值，下面我们就通过Matlab计算的模值跟arm_cfft_f32计算的模值做对比。
对比前需要先将串口打印出的数据加载到Matlab中，并给这个数组起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：
Fs = 1000; % 采样率
N= 1024; % 采样点数
n= 0:N-1; % 采样序列
t= 0:1/Fs:1-1/Fs; % 时间序列
f = n * Fs / N; %真实的频率
x = sin(2*pi*50*t) ;%原始信号
y = fft(x, N); %对原始信号做FFT变换
subplot(2,1,2);
plot(f, abs(y)); %绘制幅频相应曲线
title('Matlab计算结果');
xlabel('频率');
ylabel('幅度');
subplot(2,1,1);
plot(f,sampledata); %绘制幅频相应曲线
title('复数FFT计算结果');
xlabel('频率');
ylabel('幅度');
Matlab运行结果如下：
从上面的对比结果中可以看出，Matlab和函数arm_cfft_f32计算的结果基本是一直的。 30.2 复数FFT—基2算法
30.2.1 arm_cfft_radix2_f32
此函数已经不推荐使用，后面的版本会被删除，故不做介绍。
30.2.2 arm_cfft_radix2_q31
函数定义如下：
void arm_cfft_radix2_q31(
const arm_cfft_radix2_instance_q31 * S,
q31_t * pSrc)
参数定义：
*S points to an instance of the fixed-point CFFT/CIFFT structure.
*pSrc points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place.
注意事项：
结构const arm_cfft_radix2_instance_q31的定义如下（在文件arm_math.h文件）：
typedef struct
{
uint16_t fftLen;
uint8_t ifftFlag;
uint8_t bitReverseFlag;
q31_t *pTwiddle;
uint16_t *pBitRevTable;
uint16_t twidCoefModifier;
uint16_t bitRevFactor;
} arm_cfft_radix2_instance_q31;
下面通过在开发板上运行这个函数并计算幅频相应，然后再与Matlab计算的结果做对比。
q31_t testInput_radix2_q31_50hz;
q31_t testOutputQ31;
/*
*********************************************************************************************************
* 函 数 名: arm_cfft_radix2_q31_app
* 功能说明: 调用函数arm_cfft_radix2_q31计算幅频。
* 形 参：无
* 返 回 值: 无
*********************************************************************************************************
*/
void arm_cfft_radix2_q31_app(void)
{
uint16_t i,j;
arm_cfft_radix2_instance_q31 S;
fftSize = 1024;
ifftFlag = 0;
doBitReverse = 1;
/* 初始化结构S */
arm_cfft_radix2_init_q31(&S, fftSize, ifftFlag, doBitReverse);
/* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
for(i=0; i<1024; i++)
{
testInput_radix2_q31_50hz = 0;
/* 51.2Hz正弦波，采样率1024Hz。
arm_sin_q31输入参数的范围0-2^31， 这里每20次为一个完整的正弦波，
2^31 / 20 = 107374182.4
*/
j = i % 20;
testInput_radix2_q31_50hz = arm_sin_q31(107374182*j);
printf("%drn", testInput_radix2_q31_50hz);
}
/* 输出结果分割线 */
printf("**************************************************rn");
printf("**************************************************rn");
/* 计算CFFT */
arm_cfft_radix2_q31(&S, testInput_radix2_q31_50hz);
/* 计算模值 */
arm_cmplx_mag_q31(testInput_radix2_q31_50hz, testOutputQ31, fftSize);
/* 串口打印求解的模值 */
for(i=0; i<1024; i++)
{
printf("%drn", testOutputQ31);
}
}
运行如上函数可以通过串口打印出原始信号和计算的模值，下面我们就通过Matlab计算的模值跟arm_cfft_radix2_q31计算的模值做对比。
对比前需要先将串口打印出的数据加载到Matlab中，原始信号起名signal，函数arm_cmplx_mag_q31计算的模值起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：
Fs = 1024; % 采样率
N= 1024; % 采样点数
n= 0:N-1; % 采样序列
f = n * Fs / N; %真实的频率
y = fft(signal, N); %对原始信号做FFT变换
subplot(2,1,2);
plot(f, abs(y)); %绘制幅频相应曲线
title('Matlab计算结果');
xlabel('频率');
ylabel('幅度');
subplot(2,1,1);
plot(f,sampledata); %绘制幅频相应曲线
title('复数FFT计算结果');
xlabel('频率');
ylabel('幅度');
Matlab运行结果如下：
从上面的对比结果中可以看出，Matlab和函数arm_cfft_radix2_q31计算的频率点基本是一致的，而幅值大小不一样是因为调用函数arm_cmplx_mag_q31对数据结果做了移位处理。
30.2.3 arm_cfft_radix2_q15
函数定义：
void arm_cfft_radix2_q15(
const arm_cfft_radix2_instance_q15 * S,
q15_t * pSrc)
参数定义：
*S points to an instance of the fixed-point CFFT/CIFFT structure.
*pSrc points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place.
注意事项：
结构const arm_cfft_radix2_instance_q15的定义如下（在文件arm_math.h文件）：
typedef struct
{
uint16_t fftLen;
uint8_t ifftFlag;
uint8_t bitReverseFlag;
q15_t *pTwiddle;
uint16_t *pBitRevTable;
uint16_t twidCoefModifier;
uint16_t bitRevFactor;
} arm_cfft_radix2_instance_q15;
下面通过在开发板上运行这个函数并计算幅频相应，然后再与Matlab计算的结果做对比。
q15_t testInput_radix2_q15_50hz;
q15_t testOutputQ15;
/*
*********************************************************************************************************
* 函 数 名: arm_cfft_radix2_q15_app
* 功能说明: 调用函数arm_cfft_radix2_q15计算幅频。
* 形 参：无
* 返 回 值: 无
*********************************************************************************************************
*/
void arm_cfft_radix2_q15_app(void)
{
uint16_t i,j;
arm_cfft_radix2_instance_q15 S;
fftSize = 1024;
ifftFlag = 0;
doBitReverse = 1;
/* 初始化结构S */
arm_cfft_radix2_init_q15(&S, fftSize, ifftFlag, doBitReverse);
/* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
for(i=0; i<1024; i++)
{
/* 虚部全部置0 */
testInput_radix2_q15_50hz = 0;
/* 51.2Hz正弦波，采样率1024Hz。
arm_sin_q15输入参数的范围[0, 32768)， 这里每20次为一个完整的正弦波，
32768 / 20 = 1638.4
*/
j = i % 20;
testInput_radix2_q15_50hz = arm_sin_q15(1638*j);
printf("%drn", testInput_radix2_q15_50hz);
}
/* 输出结果分割线 */
printf("**************************************************rn");
printf("**************************************************rn");
/* 计算CFFT */
arm_cfft_radix2_q15(&S, testInput_radix2_q15_50hz);
/* 计算模值 */
arm_cmplx_mag_q15(testInput_radix2_q15_50hz, testOutputQ15, fftSize);
/* 串口打印求解的模值 */
for(i=0; i<1024; i++)
{
printf("%drn", testOutputQ15);
}
}
运行如上函数可以通过串口打印出原始信号和计算的模值，下面我们就通过Matlab计算的模值跟arm_cfft_radix2_q15计算的模值做对比。
对比前需要先将串口打印出的数据加载到Matlab中，原始信号起名signal，函数arm_cmplx_mag_q15计算的模值起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：
Fs = 1024; % 采样率
N= 1024; % 采样点数
n= 0:N-1; % 采样序列
f = n * Fs / N; %真实的频率
y = fft(signal, N); %对原始信号做FFT变换
subplot(2,1,2);
plot(f, abs(y)); %绘制幅频相应曲线
title('Matlab计算结果');
xlabel('频率');
ylabel('幅度');
subplot(2,1,1);
plot(f,sampledata); %绘制幅频相应曲线
title('复数FFT计算结果');
xlabel('频率');
ylabel('幅度');
Matlab运行结果如下：
从上面的对比结果中可以看出，Matlab和函数arm_cfft_radix2_q15计算的频率点基本是一致的，而幅值大小不一样是因为调用函数arm_cmplx_mag_q15对数据结果做了移位处理。 30.3 复数FFT—基4算法
如果希望直接调用FFT程序计算IFFT，可以用下面的方法：
30.3.1 arm_cfft_radix4_f32
此函数已经不推荐使用，后面的版本会被删除，故不做介绍。
30.3.2 arm_cfft_radix4_q31
函数定义如下：
void arm_cfft_radix4_q31(
const arm_cfft_radix4_instance_q31 * S,
q31_t * pSrc)
参数定义：
*S points to an instance of the fixed-point CFFT/CIFFT structure.
*pSrc points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place.
注意事项：
1. 结构const arm_cfft_radix4_instance_q31的定义如下（在文件arm_math.h文件）：
typedef struct
{
uint16_t fftLen;
uint8_t ifftFlag;
uint8_t bitReverseFlag;
q31_t *pTwiddle;
uint16_t *pBitRevTable;
uint16_t twidCoefModifier;
uint16_t bitRevFactor;
} arm_cfft_radix4_instance_q31;
2. 为了防止数据饱和，每次蝶形运行的结果都要除以2，故不同的长度的FFT运算的最终结果输出格式不同。具体信息如下：
Q31 CFFT
Q31 CIFFT
下面通过在开发板上运行这个函数并计算幅频相应，然后再与Matlab计算的结果做对比。
q31_t testInput_radix4_q31_50hz;
/*
*********************************************************************************************************
* 函 数 名: arm_cfft_radix4_q31_app
* 功能说明: 调用函数arm_cfft_radix4_q31_app计算幅频。
* 形 参：无
* 返 回 值: 无
*********************************************************************************************************
*/
void arm_cfft_radix4_q31_app(void)
{
uint16_t i,j;
arm_cfft_radix4_instance_q31 S;
fftSize = 1024;
ifftFlag = 0;
doBitReverse = 1;
/* 初始化结构S */
arm_cfft_radix4_init_q31(&S, fftSize, ifftFlag, doBitReverse);
/* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
for(i=0; i<1024; i++)
{
testInput_radix4_q31_50hz = 0;
/* 51.2Hz正弦波，采样率1024Hz。
arm_sin_q31输入参数的范围0-2^31， 这里每20次为一个完整的正弦波，
2^31 / 20 = 107374182.4
*/
j = i % 20;
testInput_radix4_q31_50hz = arm_sin_q31(107374182*j);
printf("%drn", testInput_radix4_q31_50hz);
}
/* 输出结果分割线 */
printf("**************************************************rn");
printf("**************************************************rn");
/* 计算CFFT */
arm_cfft_radix4_q31(&S, testInput_radix4_q31_50hz);
/* 计算模值 */
arm_cmplx_mag_q31(testInput_radix4_q31_50hz, testOutputQ31, fftSize);
/* 串口打印求解的模值 */
for(i=0; i<1024; i++)
{
printf("%drn", testOutputQ31);
}
}
运行如上函数可以通过串口打印出原始信号和计算的模值，下面我们就通过Matlab计算的模值跟arm_cfft_radix4_q31计算的模值做对比。
对比前需要先将串口打印出的数据加载到Matlab中，原始信号起名signal，函数arm_cmplx_mag_q31计算的模值起名sampledata，加载
方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：
Fs = 1024; % 采样率
N= 1024; % 采样点数
n= 0:N-1; % 采样序列
f = n * Fs / N; %真实的频率
y = fft(signal, N); %对原始信号做FFT变换
subplot(2,1,2);
plot(f, abs(y)); %绘制幅频相应曲线
title('Matlab计算结果');
xlabel('频率');
ylabel('幅度');
subplot(2,1,1);
plot(f,sampledata); %绘制幅频相应曲线
title('复数FFT计算结果');
xlabel('频率');
ylabel('幅度');
Matlab运行结果如下：
从上面的对比结果中可以看出，Matlb和函数arm_cfft_radix4_q31计算的频率点基本是一致的，而幅值大小不一样是因为调用
函数arm_cmplx_mag_q31和arm_cfft_radix4_q31对数据结果做了移位处理。
30.3.3 arm_cfft_radix4_q15
函数定义：
void arm_cfft_radix4_q15(
const arm_cfft_radix4_instance_q15 * S,
q15_t * pSrc)
参数定义：
*S points to an instance of the fixed-point CFFT/CIFFT structure.
*pSrc points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place.
注意事项：
1. 结构const arm_cfft_radix4_instance_q15的定义如下（在文件arm_math.h文件）：
typedef struct
{
uint16_t fftLen;
uint8_t ifftFlag;
uint8_t bitReverseFlag;
q15_t *pTwiddle;
uint16_t *pBitRevTable;
uint16_t twidCoefModifier;
uint16_t bitRevFactor;
} arm_cfft_radix4_instance_q15;
2. 为了防止数据饱和，每次蝶形运行的结果都要除以2，故不同的长度的FFT运算的最终结果输出格式不同。具体信息如下：
Q15 CFFT
Q15 CIFFT
下面通过在开发板上运行这个函数并计算幅频相应，然后再与Matlab计算的结果做对比。
q15_t testInput_radix4_q15_50hz;
/*
*********************************************************************************************************
* 函 数 名: arm_cfft_radix4_q15_app
* 功能说明: 调用函数arm_cfft_radix4_q15计算幅频。
* 形 参：无
* 返 回 值: 无
*********************************************************************************************************
*/
void arm_cfft_radix4_q15_app(void)
{
uint16_t i,j;
arm_cfft_radix4_instance_q15 S;
fftSize = 1024;
ifftFlag = 0;
doBitReverse = 1;
/* 初始化结构S */
arm_cfft_radix4_init_q15(&S, fftSize, ifftFlag, doBitReverse);
/* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
for(i=0; i<1024; i++)
{
/* 虚部全部置0 */
testInput_radix4_q15_50hz = 0;
/* 51.2Hz正弦波，采样率1024Hz。
arm_sin_q15输入参数的范围[0, 32768)， 这里每20次为一个完整的正弦波，
32768 / 20 = 1638.4
*/
j = i % 20;
testInput_radix4_q15_50hz = arm_sin_q15(1638*j);
printf("%drn", testInput_radix4_q15_50hz);
}
/* 输出结果分割线 */
printf("**************************************************rn");
printf("**************************************************rn");
/* 计算CFFT */
arm_cfft_radix4_q15(&S, testInput_radix4_q15_50hz);
/* 计算模值 */
arm_cmplx_mag_q15(testInput_radix4_q15_50hz, testOutputQ15, fftSize);
/* 串口打印求解的模值 */
for(i=0; i<1024; i++)
{
printf("%drn", testOutputQ15);
}
}
运行如上函数可以通过串口打印出原始信号和计算的模值，下面我们就通过Matlab计算的模值跟arm_cfft_radix4_q15计算的模值做对比。
对比前需要先将串口打印出的数据加载到Matlab中，原始信号起名signal，函数arm_cmplx_mag_q15计算的模值起名sampledata，加载
方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：
Fs = 1024; % 采样率
N= 1024; % 采样点数
n= 0:N-1; % 采样序列
f = n * Fs / N; %真实的频率
y = fft(signal, N); %对原始信号做FFT变换
subplot(2,1,2);
plot(f, abs(y)); %绘制幅频相应曲线
title('Matlab计算结果');
xlabel('频率');
ylabel('幅度');
subplot(2,1,1);
plot(f,sampledata); %绘制幅频相应曲线
title('复数FFT计算结果');
xlabel('频率');
ylabel('幅度');
Matlab运行结果如下：
从上面的对比结果中可以看出，Matlab和函数arm_cfft_radix4_q15计算的频率点基本是一致的，而幅值大小不一样是因为调用函数arm_cmplx_mag_q15
和arm_cfft_radix4_q15对数据结果做了移位处理。 30.4总结
本章节内容比较多，有兴趣的可以深入了解源码的实现。
页:
[1]