当前位置： 首页 > article >正文

MATLAB、FPGA、STM32中调用FFT计算频率、幅值及相位差

article 2024/11/14 21:27:50

系列文章目录

文章目录

  • 系列文章目录
  • 前言
  • MATLAB
  • STM32
    • 调用DSP
    • STM32中实现FFT
    • 关于初相位
  • FPGA

前言

最近在学习如何在STM32中调用FFT

MATLAB

首先对FFT进行一下说明，我们输入N个点的数据到FFT中，FFT会返回N个点的数据，这些数据都是复数，其模值就是我们用来计算频率和幅值，模值越大代表该频率占比越多，模值/N*2就是幅值。每个复数的角度代表了该频率的相位差（这里的初相位不一定准确，只有特定条件下是初相位）。FFT的结果一般是中心对称的，只需要看左半部分就行。

例如：

t = 0:1:255;
x = 30*cos(2*pi/10*t)
figure
plot(x);
y = fft(x)
figure
plot(abs(y))
tol = 2e+3;  //将不想要的频率的相位筛掉
y(abs(y) < tol) = 0;
theta = angle(y);
figure
plot(theta/pi*180);

我们看到第27个点的幅值是最大的，那么该点该表的频率就是我们所求的，假设采样率为X，那么该点的频率为X27/采样点数，该点的幅值=2930/2562，该点的相位这里显示为-71度，但很明显我们的初相位是0
在这里插入图片描述

这里我们不变采样点数，移相30度

x = 30*cos(2*pi/10*t-30/180*pi)

频率没变
在这里插入图片描述
可以看到相位确实变了30度，这就是为什么说计算不了初相位，但是对于同一频率，同样采样点数下的两个波，可以计算相位差的原因
在这里插入图片描述
那么什么情况下计算的结果是初相位呢，只有FFT的输入数据刚好是整数个周期的时候，得到的才是初相位25

x = 30*cos(2*pi/25.6*t-30/180*pi)

这里输入256个数据，那么我们设计为刚好输入10个周期的情况，很明显，得到的-30度就是我们的初相位
在这里插入图片描述

STM32

调用DSP

首先添加DSP的环境
在这里插入图片描述
之后，添加对应的.lib文件
在这里插入图片描述
添加需要的头文件
在这里插入图片描述
添加CMSIS\DSP\include的路径
在这里插入图片描述
之后就可以顺利编译了

STM32中实现FFT

#include "arm_math.h"
#include "arm_const_structs.h"

// FFT 相关参数的定义
#define FFT_SIZE 256
#define SAMPLING_FREQUENCY 10000000
float32_t inputSignal[FFT_SIZE*2];// FFT 输入信号数组
float32_t fftOutput[FFT_SIZE];// FFT 输出数组
uint32_t index_;// 存放 FFT 输出中最大值的索引
float32_t maxValue;
float32_t Vpp;
float32_t frequency ;// 用于存放计算结果的频率变量
float phase;//角度
#define     M_PI            (3.1415926f)

void fftCalculate(void)// FFT 计算函数
{
	uint8_t i;
	
    arm_cfft_f32(&arm_cfft_sR_f32_len256, inputSignal, 0, 1);// 执行 FFT 计算//这个就是快速傅里叶变换的主要接口，第一个参数可以理解为你输入到FFT里的采样点的个数；第二个参数为输入数组；第三个参数为正反变换，一般使用填0；第四个参数为位反转使能，一般使用填1。输出会覆盖掉输入
    arm_cmplx_mag_f32(inputSignal, fftOutput, FFT_SIZE );// 计算 FFT 输出的幅度，输入inputSignal，输出fftOutput
    index_ = 0;// 查找 FFT 输出中的最大值
    maxValue =  fftOutput[1]; //跳过直流分量	
	for (uint32_t i = 1; i < FFT_SIZE /2; i++)  //这里的数组第0个数据放的是直流分量，跳过就行
    {
         if (fftOutput[i] > maxValue)
        {        
           maxValue = fftOutput[i];
            index_ = i;
        }
    }
	phase = atan2(inputSignal[2*index_+1],inputSignal[2*index_]) * 180 / 3.1415926f;	//相位
    frequency = (float32_t)index_ * (float32_t)SAMPLING_FREQUENCY / (float32_t)FFT_SIZE;// 根据最大值的索引计算信号的频率
	Vpp = maxValue * 2 / 256; //幅值
}

第一个参数可以理解为你输入到FFT里的采样点个数，除了256还有512等等；第二个参数为输入数组；第三个参数为正反变换，一般使用填0；第四个参数为位反转使能，一般使用填1。输出会覆盖掉输入，所以最后计算出的复数会写入到inputSignal中。后面计算相位就需要用到这个数据

arm_cfft_f32(&arm_cfft_sR_f32_len256, inputSignal, 0, 1);

首先要注意，输入到arm_cfft_f32中的数据inputSignal是由一个实部一个虚部组成的，所以调用的时候要手动将虚部设置为0

for(i=0;i<FFT_SIZE;i++)
{
	inputSignal[i*2] = (float)SPI_AD_DATA.RX_DATA_BUF[i+503] - 128; //2+500+2
	inputSignal[i*2+1] = 0;
}

最后对计算出的相位要进行处理，因为FFT如果选择的数据是非整数个周期，计算出的相位是有一定偏移，但是在同一频率下，同一FFT计算点数下，这个偏移量是固定的，所以可以计算两个之间的相位差，但无法得到相位。

相位差就把得到的两个相位做差就行，但需要把负相位加360度转到正，最后计算出的相位差还要取绝对值

关于初相位

前面MATLAB介绍了初相位的问题，那么STM32中不能计算大部分频率的相位吗，其实有一种方法，如果能改变ADC的采样率，因为频率计算不会出现初相位这样问题，那么得到频率后，更改ADC的采样率，使得在我们设计的采样点数下，采集整数个周期，就能得到初相。

但感觉实现起来有难度

FPGA

通过调用Quaruts的官方IP核就能实现FFT名单时需要写一个外部控制模块，这里我利用ROM来进行测试，可以取intel官网找他们的IP手册
在这里插入图片描述在这里插入图片描述

在这里插入图片描述
Length：FFT变换长度64、128、256、512、1024、2048、4096、8192、16384、32768或65536。可变流还允许8、16、32、131072和262144。

Direction：傅里叶变换的方向，选择Forward：FFT（快速傅里叶变换），Reverse：FFT（快速傅里叶反变换），Bi-directional（用户可通过输入控制是FFT还是IFFT）

Calculation：计算的延迟
Throughput Latency：处理延迟

Burst：突发架构，需要的内存资源最小，平均吞吐量最低。
Buffered Brust：缓冲突发架构需要的内存资源比流式低，平均吞吐量较低。
Streaming：流式架构可以连续变换处理。
Variable Streaming：可变流式架构可以连续处理，并且可以运行时控制变换长度。在线改变FFT的大小，速度和流式差不多。前三种模式运算速度依次增大，占用资源也一次增大。
Input Otder：输入数据的顺序（顺序模式或逆序模式）
Output Order：输出数据的顺序（顺序模式或逆序模式）

Representation：数据结构和旋转因子
Block Floating（块浮点）：应用Burst、Buffered Brust、Streaming
Fixed Point(定点)和Single Floating Point（单浮点）：用于Variable Streaming
块浮点就是在数据的一帧数据中有一个共同的缩放因子，当一帧数据中有大有小的时候，共用一个缩放因子会造成误差增大。
Data Input Width：输入数据的数据宽度，8, 10, 12, 14, 16, 18, 20,24, 28, 32。
Twiddle Width：旋转因子的数据宽度，旋转因子的数据宽度不能大于输入数据的数据宽度。
Data Output Width:：输出数据的数据宽度。FFT的计算结果是输出的实部和虚部与缩放因子(EXP)的结合，缩放因子为负，输出数据需要左移（增大），为正则右移，输出的实部和虚部，缩放因子都是有符号数。
在这里插入图片描述
clk : (in) 时钟信号
reset_n : (in)复位信号，低电平有效，至少一个时钟周期
inverser : (in)低电平为FFT，高电平为IFFT
sink_valid : (in)输入数据有效信号
sink_sop : (in)输入数据起始信号，与第一个数据对齐，只需保持一个时钟周期即可
sink_eop : (in)输入数据结束信号，与最后一个数据对齐，只需保持一个时钟周期
sink_real : (in)输入数据的实部
sink_imag : (in)输入数据的虚部，一般直接置0
sink_ready : (out) IP核准备好接收数据了，实测一直高电平
sink_error : (in)输入错误信号，置0即可，不会影响。00 = no error，01 = missing start of packet (SOP)，10 = missing end of packet (EOP)，11 = unexpected EOP

source_error : (out) 输出错误信号,若输入的数据格式有误，则不进行FFT变换，并给出错误值。
source_ready : (in)输入数据准备好，置1即可，随时可以输出数据
source_sop : (out)输出数据起始信号，与输出的第一个数据对齐
source_eop : (out)输出数据的终止信号，与输出的最后一个数据对齐
source_real : (out)输出数据的实部
source_imag : (out)输出数据的虚部
source_exp : (out)数据数据的缩放因子，只有流式、突发和缓冲突发模式有。
source_vaild : (out) IFFT控制线，FFT完成时，信号置高，输出数据

输入数据时序：
在这里插入图片描述
输出数据时序：
在这里插入图片描述
多次调用时序：
在这里插入图片描述

module FFT_control
(

	input                         clk,
    input                         rst_n,
    input   [15:0]          source_real,
    input   [15:0]          source_imag,
    input   [5:0]            source_exp,   
    input                        source_valid,            
    input    [1:0]           source_error, 
    input                        source_sop, 
    input                        source_eop, 
    
    input                         sink_ready,
    
    output   reg            sink_valid,
    output   reg            sink_sop,
    output   reg            sink_eop,
    output   wire    [1:0]         sink_error,
    output   wire   [15:0]      sink_imag,                      //实部数据采集模块输入，虚部直接置0

    
    output   wire          source_ready,
    
    output   wire            inverse,    
    output   reg        	[ 10:0] 	ROM_address,
    output   wire        	[ 11:0] 	FFT_out,
    output   wire        	[ 15:0] 	sqrt_data_out
    
);

assign FFT_out=sqrt_data_out;
localparam frames_FFT=13'd1024-1'd1;
/*输入数据帧数:1024或2048看IP核内部使用*/ 
reg  [12:0] frames_in;

reg  FFT_output_start;

reg 		[ 15:0] 	fft_real_out;
reg 		[ 15:0] 	fft_image_out;
reg 		[ 31:0] 	fft_data_out;
reg 		[ 5:0] 	    fft_exp_out;

reg   FFT_count_output;
reg   FFT_count_output_MAX;
reg   FFT_count_output_start;

/* FFT输入数据虚部 */
assign  sink_imag=8'd0;
/* FFT控制信号，1:FFT，0:FFT */
assign 	inverse = 1'b0;

/* 输入错误信号 */
assign 	sink_error = 2'b0;
/*置1即可，FFT随时可以输出数据*/
assign 	source_ready=1'd1;

reg [2:0]state;
localparam 	start_state = 3'b00;
localparam 	 start_next_state= 3'b01;	
localparam 	 end_input_state= 3'b10;
localparam 	 wait_output_state= 3'b11;		


/*数据读取部分,FFT数据输入*/
always @(posedge clk or negedge rst_n) begin    
  	if(rst_n== 1'd0)begin									
		    sink_valid<=1'd0;
            sink_sop<=1'd0;
            sink_eop <= 1'd0;
            state<=start_state;
	end
	
	else begin
            case (state)
                start_state:begin		
                        sink_valid<=1'd1;           
                        sink_sop<=1'd1;
                        state<=start_next_state;
                end
                    
                start_next_state:begin 
                        sink_sop<=1'd0;
                         if(frames_in == frames_FFT-1'd1) begin
                                state <=end_input_state;		
                         end
                end
                end_input_state:begin
                         sink_valid<=1'd0;   
                         sink_eop <= 1'd1;  
                          state <=wait_output_state;
                end
                 wait_output_state:begin 
                         sink_eop <= 1'd0;  
                         FFT_output_start=1'd1;
                end
            endcase
    end
end

/*计算输入数据帧数*/
always @ (posedge clk or negedge rst_n)  begin

       if(rst_n== 1'd0) begin
             frames_in <= 'd0;
         end
         
        else begin
            if(sink_valid==1'd1) begin              //这里sink_ready会自己拉低
                    frames_in <= frames_in+1'd1;
               end
             else begin
                   frames_in <= frames_in;
               end
        end                
end

/*FFT输出数据读取*/

/* FFT输出的实部信号,数据是有符号数据，需要转换 */
always @ (*)
begin
	if(source_valid==1'd1)
		fft_real_out = source_real[15] ? (~source_real[15:0]+1) : source_real;
	else
		fft_real_out = 0; 
end

/* 相当组合逻辑，FFT输出的虚部信号 */
always @ (*)
begin
	if(source_valid==1'd1)
		fft_image_out = source_imag[15] ? (~source_imag[15:0]+1) : source_imag;
	else
		fft_image_out = 0;
end

/* 相当组合逻辑，FFT输出的数据转换 ,这个数据不能直接用来就平方根，老不稳定，很奇怪，必须加时序，加了就会延后1拍*/
always @ (posedge clk )
begin
	if(source_valid==1'd1)
        fft_data_out <= fft_real_out*fft_real_out + fft_image_out*fft_image_out;
	else
		fft_data_out <= 0;
end


/*相当组合逻辑，FFT输出的指数信号,旋转因子 */
always @ (*)
begin
	if(source_valid==1'd1)
		fft_exp_out = source_exp;
	else
		fft_exp_out = fft_exp_out;
end

/* 平方根模块IP核,fft_data_out时序电路延迟1拍 */
ALtearsqrt      SQRT_Module		 
(
	.radical 	(fft_data_out ),
	.q 			(sqrt_data_out),
);

always @(posedge clk or negedge rst_n) begin    
       if(rst_n== 1'd0) begin
             ROM_address <= 'd0;
         end
         
        else begin
                if(ROM_address<11'd1022)
                     ROM_address<=ROM_address+1'd1;
                else
                      ROM_address<='d0;
        end            
end 
endmodule

之后用signaltap抓取数据进行测试：

以Burst模式、1024数据长度、16Bits数据宽度、ROM为输入的模式进行测试（其他数据宽度经过测试效果不佳）
在这里插入图片描述
此频率代指一段FFT数据（1024）中有几个周期正弦波，因为存在时序打拍，所以这里峰值点减1为实际点位
在这里插入图片描述


http://www.kler.cn/a/308030.html

相关文章：

  • [ 网络安全介绍 5 ] 为什么要学习网络安全？
  • 图片画廊 day2 （可复制源码）
  • matlab建模入门指导
  • request爬虫库的小坑
  • 《情商》提升：增强自我意识，学会与情绪共处
  • 10款翻译工具实践体验感受与解析！！！！！
  • 深度学习----------------------文本预处理
  • 进阶SpringBoot之异步任务、邮件任务和定时执行任务
  • 2024最新版MySQL详细学习教程
  • 算法两道题
  • PyCharm 安装
  • 【重学 MySQL】二十八、SQL99语法新特性之自然连接和 using 连接
  • 无人机在战争方面的应用！！！
  • 计算机毕业设计 基于协同过滤算法的个性化音乐推荐系统 Java+SpringBoot+Vue 前后端分离 文档报告 代码讲解 安装调试
  • unity3d入门教程七
  • 如何编写智能合约——基于长安链的Go语言的合约开发
  • 我想要抓取新加坡当地电商平台数据，使用什么地区的IP最合适
  • C++数据排序( 附源码 )
  • linux-硬件与设备管理-硬件信息查看
  • 计算机网络 第3章 数据链路层
  • JS - 获取剪切板内容 Clipboard API
  • 数据结构——栈和队列（队列的定义、顺序队列以及链式队列的基本操作）
  • opencv学习：图像直方图均衡化与对比度受限的自适应直方图均衡化及实验代码
  • 针对特定接口记录审核日志类的写入数据库的方法
  • raksmart的G口大流量服务器怎么样?
  • C++学习， 数据抽象