【安富莱DSP教程】第32章 实数FFT的实现
特别说明：完整45期数字信号处理教程，原创高性能示波器代码全开源地址：链接第32章 实数FFT的实现
本章主要讲解实数的浮点和定点Q31，Q15的实现。关于这部分的知识点和函数的计算结果上，官方的文档有一些小错误，在章节中会跟大家详细讲述，还有一个要注意的问题，调用实数FFT函数一定要使用CMSIS-DSP V1.4.4及其以上版本，以前的版本有bug。
本章节使用的复数FFT函数来自ARM官方库的TransformFunctions部分
32.1 复数FFT
32.2 复数FFT-基2算法
32.3 复数FFT-基4算法
32.4 总结
32.1 实数FFT
32.1.1 描述
CMSIS DSP库里面包含一个专门用于计算实数序列的FFT库，很多情况下，用户只需要计算实数序列即可。计算同样点数FFT的实数序列要比计算同样点数的虚数序列有速度上的优势。
快速的rfft算法是基于混合基cfft算法实现的。
一个N点的实数序列FFT正变换采用下面的步骤实现：
由上面的框图可以看出，实数序列的FFT是先计算N/2个实数的CFFT，然后再重塑数据进行处理从而获得半个FFT频谱即可（利用了FFT变换后频谱的对称性）。
一个N点的实数序列FFT逆变换采用下面的步骤实现：
实数FFT支持浮点，Q31和Q15三种数据类型。
32.2 实数FFT
32.2.1 arm_rfft_fast_f32
函数定义如下：
void arm_rfft_fast_f32(
arm_rfft_fast_instance_f32 * S,
float32_t * p, float32_t * pOut,
uint8_t ifftFlag)
参数定义：
*S points to an arm_rfft_fast_instance_f32 structure.
*p points to the input buffer.
*pOut points to the output buffer.
ifftFlag RFFT if flag is 0, RIFFT if flag is 1
注意事项：
结构arm_rfft_fast_instance_f32的定义如下（在文件arm_math.h文件）：
typedef struct
{
arm_cfft_instance_f32 Sint; /**< Internal CFFT structure. */
uint16_t fftLenRFFT; /**< length of the real sequence */
float32_t * pTwiddleRFFT; /**< Twiddle factors real stage*/
} arm_rfft_fast_instance_f32 ;
下面通过在开发板上运行函数arm_rfft_fast_f32和arm_cfft_f32计算幅频响应，然后将相应的频率响应结果在Matlab上面绘制出来。
/*
*********************************************************************************************************
* 函 数 名: arm_rfft_fast_f32_app
* 功能说明: 调用函数arm_rfft_fast_f32计算1024点实数序列的幅频响应并跟使用函数arm_cfft_f32计算结果做对比
* 形 参：无
* 返 回 值: 无
*********************************************************************************************************
*/
static void arm_rfft_fast_f32_app(void)
{
uint16_t i;
arm_rfft_fast_instance_f32 S;
/* 实数序列FFT长度 */
fftSize = 1024;
/* 正变换 */
ifftFlag = 0;
/* 初始化结构体S中的参数 */
arm_rfft_fast_init_f32(&S, fftSize);
/* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
for(i=0; i<1024; i++)
{
/* 50Hz正弦波，采样率1KHz */
testInput_f32_10khz = 1.2f*arm_sin_f32(2*3.1415926f*50*i/1000)+1;
}
/* 1024点实序列快速FFT */
arm_rfft_fast_f32(&S, testInput_f32_10khz, testOutput_f32_10khz, ifftFlag);
/* 为了方便跟函数arm_cfft_f32计算的结果做对比，这里求解了1024组模值，实际函数arm_rfft_fast_f32
只求解出了512组
*/
arm_cmplx_mag_f32(testOutput_f32_10khz, testOutput, fftSize);
/* 串口打印求解的模值 */
for(i=0; i<fftSize; i++)
{
printf("%f\r\n", testOutput);
}
printf("****************************分割线***************************************\r\n");
for(i=0; i<1024; i++)
{
/* 虚部全部置零 */
testInput_f32_10khz = 0;
/* 50Hz正弦波，采样率1KHz ,作为实部 */
testInput_f32_10khz = 1.2f*arm_sin_f32(2*3.1415926f*50*i/1000)+1;
}
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<fftSize; i++)
{
printf("%f\r\n", testOutput);
}
}
运行如上函数可以通过串口打印出函数arm_rfft_fast_f32和arm_cfft_f32计算的幅频模值，下面通过Matlab绘制波形来对比这两种模值。
对比前需要先将串口打印出的两组数据加载到Matlab中，arm_rfft_fast_f32的计算结果起名signal，arm_cfft_f32的计算结果起名sampledata，
加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：
Fs = 1000; % 采样率
N= 1024; % 采样点数
n= 0:N-1; % 采样序列
f = n * Fs / N; %真实的频率
subplot(3,1,1);
plot(f,signal); %绘制RFFT结果
title('实数FFT');
xlabel('时间');
ylabel('幅值');
subplot(3,1,2);
plot(f,sampledata); %CFFT结果
title('复数FFT');
xlabel('时间');
ylabel('幅值');
Matlab运行结果如下：
从上面的前512点对比中，我们可以看出两者的计算结果是相符的。这里有一点要特别注意，官方文档中对于函数arm_rfft_fast_f32输出结果的实部，虚部排列顺序说明是错误的。函数arm_rfft_fast_f32的输出结果仍然是实部，虚部，实部，虚部….. 依次排列下去。
函数arm_rfft_fast_f32在计算直流分量（也就是频率为0的值）的虚部上是有错误的。关于这点大家可以将实际的实部和虚部输出结果打印出来做对比，但差别很小，基本可以忽略。
32.2.2 arm_rfft_q15
函数定义如下：
void arm_rfft_q15(
const arm_rfft_instance_q15 * S,
q15_t * pSrc,
q15_t * pDst)
参数定义：
*S points to an instance of the Q15 RFFT/RIFFT structure.
[in]*pSrc points to the input buffer.
*pDst points to the output buffer.
return none.
注意事项：
结构arm_rfft_instance_q15的定义如下（在文件arm_math.h文件）：
typedef struct
{
uint32_t fftLenReal;
uint8_t ifftFlagR;
uint8_t bitReverseFlagR;
uint32_t twidCoefRModifier;
q15_t *pTwiddleAReal;
q15_t *pTwiddleBReal;
const arm_cfft_instance_q15 *pCfft;
} arm_rfft_instance_q15;
下面通过在开发板上运行函数arm_rfft_q15和arm_cfft_f32计算幅频响应，然后将相应的频率响应结果在Matlab上面绘制出来。
/*
*********************************************************************************************************
* 函 数 名: arm_rfft_q15_app
* 功能说明: 调用函数arm_rfft_q15计算1024点实数序列的幅频响应并跟使用函数arm_cfft_f32计算的结果做对比。
* 形 参：无
* 返 回 值: 无
*********************************************************************************************************
*/
static void arm_rfft_q15_app(void)
{
uint16_t i,j;
arm_rfft_instance_q15 S;
/* 实数序列FFT长度 */
fftSize = 1024;
/* 正变换 */
ifftFlag = 0;
/* 码位倒序 */
doBitReverse = 1;
/* 初始化结构体S */
arm_rfft_init_q15(&S, fftSize, ifftFlag, doBitReverse);
/* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
for(i=0; i<1024; i++)
{
/* 51.2Hz正弦波，采样率1024Hz。
arm_sin_q15输入参数的范围[0, 32768)， 这里每20次为一个完整的正弦波，
32768 / 20 = 1638.4
*/
j = i % 20;
testInput_q15_50hz = arm_sin_q15(1638*j);
}
/* 1024点实序列快速FFT */
arm_rfft_q15(&S, testInput_q15_50hz, testOutput_q15_50hz);
/* 由于输出结果的格式是Q5，所以这里将定点数转换为浮点数 */
for(i = 0; i < fftSize; i++)
{
testOutput_f32_10khz = (float32_t)testOutput_q15_50hz/32;
}
/* 为了方便对比，这里求解了1024组复数，实际上面的变化只有512组
实际函数arm_rfft_q15只求解出了512组*/
arm_cmplx_mag_f32(testOutput_f32_10khz, testOutput, fftSize);
/* 串口打印求解的模值 */
for(i=0; i<fftSize; i++)
{
printf("%f\r\n", testOutput);
}
printf("****************************分割线***************************************\r\n");
for(i=0; i<1024; i++)
{
/* 51.2Hz正弦波，采样率1024Hz。
arm_sin_q15输入参数的范围[0, 32768)， 这里每20次为一个完整的正弦波，
32768 / 20 = 1638.4
*/
j = i % 20;
testInput_f32_10khz = (float32_t) arm_sin_q15(1638*j)/32768;
/* 虚部全部置零 */
testInput_f32_10khz = 0;
}
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<fftSize; i++)
{
printf("%f\r\n", testOutput);
}
}
运行如上函数可以通过串口打印出函数arm_rfft_q15和arm_cfft_f32计算的幅频模值，下面通过Matlab绘制波形来对比这两种模值。
对比前需要先将串口打印出的两组数据加载到Matlab中，arm_rfft_q15的计算结果起名signal，arm_cfft_f32的计算结果起名sampledata，
加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：
Fs = 1000; % 采样率
N= 1024; % 采样点数
n= 0:N-1; % 采样序列
f = n * Fs / N; %真实的频率
subplot(3,1,1);
plot(f,signal); %绘制RFFT结果
title('实数FFT');
xlabel('时间');
ylabel('幅值');
subplot(3,1,2);
plot(f,sampledata); %CFFT结果
title('复数FFT');
xlabel('时间');
ylabel('幅值');
Matlab运行结果如下：
从上面的前512点对比中，我们可以看出两者的计算结果是相符的。这里有一点要特别注意，官方文档中对于函数arm_rfft_q31输出结果的实部，虚部排列顺序说明是错误的。函数arm_rfft_q31的输出结果仍然是实部，虚部，实部，虚部….. 依次排列下去。
32.2.3 arm_rfft_q31
函数定义如下：
void arm_rfft_q31(
const arm_rfft_instance_q31 * S,
q31_t * pSrc,
q31_t * pDst)
参数定义：
*S points to an instance of the Q31 RFFT/RIFFT structure.
[in]*pSrc points to the input buffer.
*pDst points to the output buffer.
return none.
注意事项：
结构arm_rfft_instance_q31的定义如下（在文件arm_math.h文件）：
typedef struct
{
uint32_t fftLenReal;
uint8_t ifftFlagR;
uint8_t bitReverseFlagR;
uint32_t twidCoefRModifier;
q31_t *pTwiddleAReal;
q31_t *pTwiddleBReal;
const arm_cfft_instance_q31 *pCfft;
} arm_rfft_instance_q31;
下面通过在开发板上运行函数arm_rfft_q31和arm_cfft_f32计算幅频响应，然后将相应的频率响应结果在Matlab上面绘制出来。
/*
*********************************************************************************************************
* 函 数 名: arm_rfft_q31_app
* 功能说明: 调用函数arm_rfft_q31计算1024点实数序列的幅频响应并跟使用函数arm_cfft_f32计算的结果做对比。
* 形 参：无
* 返 回 值: 无
*********************************************************************************************************
*/
static void arm_rfft_q31_app(void)
{
uint16_t i,j;
arm_rfft_instance_q31 S;
/* 实数序列FFT长度 */
fftSize = 1024;
/* 正变换 */
ifftFlag = 0;
/* 码位倒序 */
doBitReverse = 1;
/* 初始化结构体S */
arm_rfft_init_q31(&S, fftSize, ifftFlag, doBitReverse);
/* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
for(i=0; i<1024; i++)
{
/* 51.2Hz正弦波，采样率1024Hz。
arm_sin_q31输入参数的范围0-2^31， 这里每20次为一个完整的正弦波，
2^31 / 20 = 107374182.4
*/
j = i % 20;
testInput_q31_50hz = arm_sin_q31(107374182*j);
}
/* 1024点实序列快速FFT */
arm_rfft_q31(&S, testInput_q31_50hz, testOutput_q31_50hz);
/* 由于输出结果的格式是Q21，所以这里将定点数转换为浮点数 */
for(i = 0; i < fftSize; i++)
{
/* 输出的数据是11.21格式，2^21 = 4194304*/
testOutput_f32_10khz = (float32_t)testOutput_q31_50hz/2097152;
}
/* 为了方便对比，这里求解了1024组复数，实际上面的变化只有512组
实际函数arm_rfft_q31只求解出了512组*/
arm_cmplx_mag_f32(testOutput_f32_10khz, testOutput, fftSize);
/* 串口打印求解的模值 */
for(i=0; i<fftSize; i++)
{
printf("%f\r\n", testOutput);
}
printf("****************************分割线***************************************\r\n");
for(i=0; i<1024; i++)
{
/* 51.2Hz正弦波，采样率1024Hz。
arm_sin_q31输入参数的范围0-2^31， 这里每20次为一个完整的正弦波，
2^31 / 20 = 107374182.4
*/
j = i % 20;
testInput_f32_10khz = (float32_t)arm_sin_q31(107374182*j)/2147483648;
/* 虚部全部置零 */
testInput_f32_10khz = 0;
}
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<fftSize; i++)
{
printf("%f\r\n", testOutput);
}
}
运行如上函数可以通过串口打印出函数arm_rfft_q15和arm_cfft_f32计算的幅频模值，下面通过Matlab绘制波形来对比这两种模值。
对比前需要先将串口打印出的两组数据加载到Matlab中，arm_rfft_q15的计算结果起名signal，arm_cfft_f32的计算结果起名sampledata，
加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：
Fs = 1000; % 采样率
N= 1024; % 采样点数
n= 0:N-1; % 采样序列
f = n * Fs / N; %真实的频率
subplot(3,1,1);
plot(f,signal); %绘制RFFT结果
title('实数FFT');
xlabel('时间');
ylabel('幅值');
subplot(3,1,2);
plot(f,sampledata); %CFFT结果
title('复数FFT');
xlabel('时间');
ylabel('幅值');
Matlab运行结果如下：
从上面的前512点对比中，我们可以看出两者的计算结果是相符的。这里有一点要特别注意，官方文档中对于函数arm_rfft_q31输出结果的实部，虚部排列顺序说明是错误的。函数arm_rfft_q31的输出结果仍然是实部，虚部，实部，虚部….. 依次排列下去。
32.3 总结
使用实数FFT计算的时候要特别的注意本章节提到的几个错误点。 在实数fft例程中，输入序列按照例程输入来看不是按照实部，虚部，实部，虚部..... 存储的吧？
/* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
for(i=0; i<1024; i++)
{
/* 50Hz正弦波，采样率1KHz */
testInput_f32_10khz = 1.2f*arm_sin_f32(2*3.1415926f*50*i/1000)+1;
}
/ confidentlin 发表于 2018-7-16 16:15
在实数fft例程中，输入序列按照例程输入来看不是按照实部，虚部，实部，虚部..... 存储的吧？
/* 按照实部 ...
是的，注释有误，忘了修改了。 那输入是按实部，实部，实部，实部.....存储的，输出是实部，虚部，实部，虚部.....存储的是吧？ confidentlin 发表于 2018-7-16 16:20
那输入是按实部，实部，实部，实部.....存储的，输出是实部，虚部，实部，虚部.....存储的是吧？
对。 eric2013 发表于 2018-7-16 16:22
对。
请问Q15的实数IFFT算法在使用的过程中，在进行完FFT运算后，直接将FFT的输出结果 testOutput_q15_50hz作为IFFT输入，可以吗？
static void arm_rifft_q15_app(void)
{
fftSize = 1024;
/* ??? */
ifftFlag = 1;
/* ???? */
doBitReverse = 1;
uint8_t i;
arm_rfft_instance_q15 S;
arm_rfft_init_q15(&S, fftSize, ifftFlag, doBitReverse);
arm_rfft_q15(&S, testOutput_q15_50hz, testOutput1_q15_50hz);
printf("printf the ifft data.");
for(i=0; i<fftSize; i++)
{
printf("%d\r\n", testOutput1_q15_50hz);
}
} liwei0225 发表于 2021-1-12 09:18
请问Q15的实数IFFT算法在使用的过程中，在进行完FFT运算后，直接将FFT的输出结果 testOutput_q15_50hz作 ...
你这个是用arm_rfft_q15做的FFT，然后将输出数据再用arm_rfft_q15做FFT逆变换？
早期我在示波器上玩过，相位上有些错位，正弦波应该好点，你用我的方法玩下
http://www.armbbs.cn/forum.php?m ... 2408&extra=page%3D1
上面是正变化，下面是逆变换
https://img.anfulai.cn/dz/attachment/forum/pw/Fid_2/2_58_0604d89c69c9a7c.jpg?77
eric2013 发表于 2021-1-12 10:27
你这个是用arm_rfft_q15做的FFT，然后将输出数据再用arm_rfft_q15做FFT逆变换？
早期我在示波器上玩过 ...
关于 arm_rfft_q15 这个api，感觉很玄乎，自己测试发现有问题，根据官网的图片如下：
arm_rfft_q15正变换：
arm_rfft_q15 反变换：
另外参考网上唯一的一个arm_rfft_q15 正变换又反变换的例子如下：
https://stackoverflow.com/questions/46147919/cmsis-real-fft-on-8192-samples-in-q15
我把主要源码拷贝如下：
if(arm_rfft_init_q15(&S, 8192, 0, 1) != ARM_MATH_SUCCESS)
return state_error;
// perform FFT
arm_rfft_q15(&S, (q15_t*)0xC0000000, (q15_t*)0xC0400000);
// Post-shift by 12, in place (see doc)
arm_shift_q15((q15_t*)0xC0400000, 12, (q15_t*)0xC0400000, 16384);
// Init inverse FFT
if(arm_rfft_init_q15(&S, 8192, 1, 1) != ARM_MATH_SUCCESS)
return state_error;
// Perform iFFT
arm_rfft_q15(&S, (q15_t*)0xC0400000, (q15_t*)0xC0800000);
// Post shift
arm_shift_q15((q15_t*)0xC0800000, 12, (q15_t*)0xC0800000, 8192);
自己尝试发现基本上经过正变换又反变换，没法复原原始信号，查看正变换之后经过12移位16bits的数据基本很多都32767了，都饱和了。
不知道FFT 和IFFT是不是上面代码那样使用，如果是为什么信号不能复原，另外IFFT之后是不是还需要shift的操作，能否指教下，谢谢！
joinee0208 发表于 2021-3-11 11:10
关于 arm_rfft_q15 这个api，感觉很玄乎，自己测试发现有问题，根据官网的图片如下：
arm_rfft_q15正变 ...
你这个帖子的楼上回复，我用的是ST的那个汇编FFT做的逆变换，效果还正常。
ARM的这个定点FFT，我还没有测试过，后面我测试下，我们交流下。你可以先用浮点FFT测试下，那个没有定标问题，测试还方便些。 eric2013 发表于 2021-3-11 17:39
你这个帖子的楼上回复，我用的是ST的那个汇编FFT做的逆变换，效果还正常。
ARM的这个定点FFT，我还没 ...
恩，之所以用q15，是因为数据都是16bits的音频，用q15刚好不用转化，如果用float的就需要每次fft之前从float转化成q15，内存扩大一倍，还消耗计算资源，索性直接用q15，我查看了内网和外网的所有例子（除了上面我提到的stackoverflow的demo），几乎没有q15正变换再反变换回来的成功的例子(基本上正变换之后就不管了)，所以觉得很奇怪，也有可能是我理解不对，后续你有测试进展，还望指教下，谢谢了！ eric2013 发表于 2021-3-11 17:39
你这个帖子的楼上回复，我用的是ST的那个汇编FFT做的逆变换，效果还正常。
ARM的这个定点FFT，我还没 ...
斑竹有没有空啊，能否帮忙测试下啊:) joinee0208 发表于 2021-5-12 08:54
斑竹有没有空啊，能否帮忙测试下啊
有时间了试试。 请问Eric，按您的例程，做1024点FFT，采样率100HZ, 调用函数 arm_rfft_fast_init_f32(&S_FFT, 1024);
arm_rfft_fast_f32(&S_FFT, m_stFFTBuffer.FFTBufIn, m_stFFTBuffer.FFTBufOut, 0);
m_stFFTBuffer.FFTBufIn 的数组长度是1024，实际每次采样256个点的数据存入，后面的全都是补零，每采样完256个数据后执行1024点的FFT, 得到的结果 m_stFFTBuffer.FFTBufOut 使用KEIL watch窗口观察 逐渐变大，25秒钟以后就溢出了，请问 是 arm_rfft_fast_f32 必须要采样满1024个点才能处理吗 ？同样的方式使用cr4_fft_1024_stm32 函数处理就不会有问题 xml2028 发表于 2021-8-16 11:18
请问Eric，按您的例程，做1024点FFT，采样率100HZ, 调用函数 arm_rfft_fast_init_f32(&S_FFT, 1024);
arm_ ...
采集256个点后面点补0没问题。
输出数据逐渐变大直到溢出，不正常。 eric2013 发表于 2021-8-16 14:09
采集256个点后面点补0没问题。
输出数据逐渐变大直到溢出，不正常。
我用的雅特力的M4内核的AT32F403A单片机，
硬汉哥的DSP教程说：相应文件添加后还有最重要一条，要把stm32_dsp.h文件中的STM32F4的头文件， 我将stm32f4xx.h"替换成了 #include "at32f4xx.h" 这个会不会影响函数arm_rfft_fast_f32的调用 ，我编译也没有报错 现在发现调用了函数 arm_rfft_fast_f32(&S_FFT, m_stFFTBuffer.FFTBufIn, m_stFFTBuffer.FFTBufOut, 0); 后，
数组m_stFFTBuffer.FFTBufIn 内的数据发生了很大变化 xml2028 发表于 2021-8-16 15:57
现在发现调用了函数 arm_rfft_fast_f32(&S_FFT, m_stFFTBuffer.FFTBufIn, m_stFFTBuffer.FFTBufOut, 0); 后 ...
这是正常的，此变量空间也要参与运算。
使用DSP库的实数FFT函数arm_rfft_fast_f32，注意输入变量也会参与运算，存入临时数值
http://www.armbbs.cn/forum.php?mod=viewthread&tid=107963&fromuid=58
(出处: 硬汉嵌入式论坛)
请问对于不定长的采样样本FFT的实现，有哪些参考案例或者库吗？
非2的N次方长度的。
感谢。 darkness 发表于 2022-2-24 21:04
请问对于不定长的采样样本FFT的实现，有哪些参考案例或者库吗？
非2的N次方长度的。
感谢。
补零就行。 joinee0208 发表于 2021-3-11 11:10
关于 arm_rfft_q15 这个api，感觉很玄乎，自己测试发现有问题，根据官网的图片如下：
arm_rfft_q15正变 ...
哥们 你搞定了么我现在也遇到逆变换的问题了
页:
[1]