【嵌入式】从Matlab到STM32F4实现频谱分析（附源代码+手推公式）

基于STM32F407单片机DSP库，实现了AD采集、FFT频谱分析等功能，并通过串口将频率向量、频谱值、采样波形发送至PC端，无需后续处理，简单好用

文章共4,021字 · 阅读需要大约14分钟

一键AI生成摘要，助你高效阅读

问答

m0_56655321

1688人浏览 · 2023-08-17 09:10:16

m0_56655321 · 2023-08-17 09:10:16 发布

1. 概要

本文记录了自己基于STM32F4单片机DSP库实现“AD采集+FFT频谱分析”的项目。在例程基础上实现了fftshift、频率换算、串口打印等功能，并对100Hz方波信号进行了相应测试，效果良好。

2. 频谱分析原理概述

FFT（快速傅里叶变换）是DFT（离散傅里叶变换）的一种快速算法，基于采样和FFT算法可以进行模拟信号的频谱分析。本次分享不介绍FFT和DFT的基本原理，而希望直接以Matlab仿真的方式展示模拟信号频谱分析的过程。

如图，设置仿真时间为 T_max=0.1s, 采样率为 fs=1000Hz, 信号频率 fm=50Hz（周期 $T_{s}=\frac{1}{f_{m}}=0.02s$ ）在仿真时间内信号共有 $\frac{T_{max}}{T_{s}}=\frac{0.1s}{0.02s}=5$ 个周期。其频谱结果如右图所示，在±50Hz频率都有尖峰，表示在50Hz频率处有较强的频率分量。

（另外，还可以看到频谱图的幅值较大，因此在实际中需要对幅度进行“归一化”，对于周期信号而言，要将fft运算结果除以采样点数N; 对于非周期信号，要将fft运算结果除以采样频率fs, 具体可见《数字信号处理》相关教材）

%% Matlab实现频谱分析
fs = 1000;    % 采样率
fm = 50;      % 信号频率
dt = 1/fs;    % 采样时间间隔
Tmax = 0.1;   % 仿真结束时间
N = floor(Tmax / dt);           % 抽样点数
t = 0 : dt : (N - 1) * dt;      % 时间向量
y = cos(2 * pi * fm * t);       % 余弦信号

figure()
plot(t,y,'LineWidth',1.2);
xlabel('时间/s')
ylabel('幅度')
title('时域波形')

df = fs / N;     % 频率间隔
f_index = [-N/2 : (N/2 - 1)] * df;   % 频率向量
freq = fftshift(fft(y));             % 频谱

figure();
plot(f_index, abs(freq),'LineWidth',1.2);
ylim([0 , 1.2 * max(abs(freq))])
xlabel('频率/Hz')
ylabel('幅度')
title('频谱图')

3. DSP库函数实现fft运算相关函数

在DSP库中，主要提供了以下几个函数实现fft运算

arm_status arm_cfft_radix4_init_f32(arm_cfft_radix4_instance_f32 * S,uint16_t fftLen,uint8_t ifftFlag,uint8_t bitReverseFlag)
void arm_cfft_radix4_f32(const arm_cfft_radix4_instance_f32 * S,float32_t * pSrc)
void arm_cmplx_mag_f32(float32_t * pSrc,float32_t * pDst,uint32_t numSamples)

函数具体解释如下(转自[1])

/*
 * 函数1  初始化FFT运算相关参数
 *	fftLen 			用于指定 FFT长度（16/64/256/1024/4096）本章设置为1024
 *	ifftFlag 		用于指定是傅里叶变换(0)还是反傅里叶变换(1) 本章设置为0
 *	bitReverseFlag  用于设置是否按位取反 本章设置为 1
 * */
arm_status arm_cfft_radix4_init_f32(
	arm_cfft_radix4_instance_f32 * S,
	uint16_t fftLen,uint8_t ifftFlag,uint8_t bitReverseFlag)


/*
 *  函数2  执行基 4 浮点FFT运算
 * 	S结构体指针参数 先由 arm_cfft_radix4_init_f32 函数设置好 然后传入该函数的
 * 	pSrc 		  传入采集到的输入信号数据（实部+虚部形式）同时FFT变换后的数据也按顺序存放在pSrc里面pSrc 				  必须大于等于2倍fftLen长度
 * */
void arm_cfft_radix4_f32(const arm_cfft_radix4_instance_f32 * S,float32_t * pSrc)


/*
 *  函数3  计算复数模值 对 FFT 变换后的结果数据，执行取模操作
 *	pSrc 		复数输入数组（大小为 2*numSamples）指针，指向 FFT 变换后的结果
 *	pDst		输出数组（大小为 numSamples）指针，存储取模后的值
 *	numSamples  就是总共有多少个数据需要取模
 * */
void arm_cmplx_mag_f32(float32_t * pSrc,float32_t * pDst,uint32_t numSamples)

简言之，基本流程是

调用函数1初始化fft算法配置，包括FFT点数、正/反FFT变换等
调用函数2函数进行fft运算，运算结果为复数形式（实部、虚部交替放置，eg. a[0] = 1, a[1] = 5, 则表示第一个复数为1+5j）
调用函数3计算fft结果幅值，运算结果为实数形式（即a[0] = 1, a[1] = 3,则有b[0] = sqrt(10)）

但是，本文作者在实际测试中发现函数3得不到相应结果，于是自定义了一个计算幅值的函数cmplx_mag。

4. 项目整体流程

4.1 相关外设初始化

在本项目中，主要涉及以下几个外设

ADC—— 用于电压采集
定时器TIM2——用于“手动控制”ADC采样速率
定时器TIM3——用于记录采集一帧（FFT采样点数）数据需要的时长，从而计算实际采样率
串口——打印数据

4.2 基于TIM2和ADC采集一帧数据

在本项目中，考虑使用“TIM2定时器中断 + ADC手动采集”的方式控制采样速率，并通过TIM3对一帧数据的采集时间进行记录，便于后续对于真实采样率进行计算，其主要代码如下：

//TIM2中断程序
if(num <= FFT_LENGTH - 1)      //采集一帧数据
{				
	if(num == 0)	    //刚开始采集
	{
		TIM_Cmd(TIM3,DISABLE);   	//T3定时器失能
		T3_over_time = 0;			//溢出clear
		TIM_SetCounter(TIM3, 0);	//计数器clear			
		TIM_Cmd(TIM3,ENABLE);   	//T3定时器使能	
	}
	fft_inputbuf[2*num] = get_advalue(); //读取ADC采集寄存器的数值（实部）
	fft_inputbuf[2*num + 1] = 0;		 //虚部设置为0
	if(num == FFT_LENGTH - 1)		//采集结束
	{
		TIM_Cmd(TIM3,DISABLE);   	//T3定时器失能
		T3_clock = TIM_GetCounter(TIM3)|(T3_over_time<<16);  //采集FFT_LENGTH个点需要TIM3的计数值
		SAMPLE_RATE = FFT_LENGTH * 1000000 /(1.0*T3_clock);	 //计算实际采样率
    }
}
num ++

在T2中断程序中，对num值进行累加。当num = 0时表示当前为采样的第一个数据点，因此需要将TIM3计数状态清零，随后将AD采集数据写入fft_inputbuf数组的实部（虚部置0）；当num = FFT_LENGTH – 1表示采集结束，计算采样频率。

对于采样频率计算，采样一帧数据（FFT点数）所花的时间

$T_{perframe}=\frac{Clk_{TIM3}}{f_{clk}}$

其中 $T_{perframe}$ 表示一帧数据采样的时间， $Clk_{TIM3}$ 表示这段时间内TIM3的计数值， $f_{clk}$ 表示TIM3的计数频率。在本项目中，设置TIM3分频系数 $TIM_{Prescaler}$ 为84，因此有

$f_{clk}=\frac{f_{sys}/2}{TIM_{Prescaler}}=\frac{84MHz}{84}=1MHz$

其中 $f_{sys}$ 为系统时钟频率168MHz。

则有采样间隔

$dT=\frac{T_{perframe}}{FFT_{LENGTH}}$

采样频率

$f_{s}=\frac{1}{dT}=\frac{FFT_{LENGTH}}{T_{perframe}}=\frac{FFT_{LENGTH\times f_{clk}}}{Clk_{TIM3}}$

4.3 FFT和幅值计算

完成数据采集之后，需要调用DSP库函数计算FFT，主要代码如下：

arm_cfft_radix4_f32(&scfft,fft_inputbuf);	    // FFT计算（基4）
cmplx_mag(fft_inputbuf,FFT_LENGTH);        		// 计算幅值
my_fftshift(FFT_LENGTH);                  		// 移位
cal_f_index(SAMPLE_RATE,FFT_LENGTH);		    // 计算频率分量

// describe: 计算复数幅度
// *src : 输入复试数据，实部 + 虚部交替排列
// numSamples: 复数数据长度 = length(scr)/2
void cmplx_mag(float *src, unsigned int numSamples)    
{
	float real;
	float imag;
	uint16_t i;
    for ( i = 0; i < numSamples; i++) 
	{
       real = src[2*i];       // 实部
       imag = src[2*i + 1];   // 虚部
       fft_outputbuf[i] = sqrtf(real*real + imag*imag);   // 幅度
    }
}

// describe: FFTSHIFT(FFT移位)
// numSamples: 采样点数（FFT点数）
void my_fftshift(unsigned int numSamples)  
{
	uint16_t i;
    for ( i = 0; i < numSamples/2; i++)    			// 前半部分向后平移
	{
		fft_output_shift[i + numSamples/2] = fft_outputbuf[i];
    }
    for ( i = numSamples/2; i < numSamples; i++)    // 后半部分向前平移
	{
		fft_output_shift[i - numSamples/2] = fft_outputbuf[i];
    }
}

// describe: 计算频率分量
// sample_rate : 实际采样率/Hz
// numSamples: 采样点数（FFT点数）
void cal_f_index(int32_t sample_rate,int numSamples)		 
{
	int16_t i;
	float df = sample_rate / (1.0*numSamples); 		// 频率分辨率
    for ( i = -(numSamples/2); i < numSamples/2; i++)  	
	{
		f_index[i + numSamples/2] = df * i;         // 计算频率向量
    }
}

首先调用arm_cfft_radix4_f32函数进行FFT运算，结果保存在fft_inputbuf复数数组中（实部加虚部交替），

其次调用cmplx_mag函数计算fft_inputbuf复数数组的赋值，结果保存在全局变量fft_outputbuf中，

然后调用my_fftshift函数对fft_outputbuf结果进行移位，结果保存在全局变量fft_output_shift中，

最后调用cal_f_index函数计算频率分辨率，结果保存在全局变量f_index中。