matlab中fft定点运算,可用于嵌入式计算的定点FFT算法 (转载)
1
时域抽取法的定点FFT程序3.0 !
简介
本示例程序演示了如何计算定点和浮点各个长度和位数的按时间抽取FFT算法。所有程序在VC6下编译通过。
最精确的FFT计算最好采用浮点,但浮点对于嵌入式设备计算量太大,有的朋友想采用定点,但网上资料不多,且长度位数要求不一,于是本人制作了这个版本。
此程序可以实现对长度为64,128,256,512,1024,2048的数据进行定点FFT(如果需要其他长度可自行修改)以及浮点版本。
根据处理器处理能力不同,对于数据类型的长度,提供了8位版本,16位版本和32位版本,可用于单片机等不同位数的嵌入式处理器。当然,三个版本主要是给定了不同的处理精度。
此程序采用了运算级间动态缩放舍去最低位,保证了输出不溢出,且提供了尽可能高的精度。具有自动适应输入值大小的能力。因此有较高的信噪比。
使用说明
FFT计算内容整合在timefft.h中,可将其加入到图形环境中或者控制台中,正如示例所示。
有两个值需要根据环境修改:
长度通过#define N 长度值
决定,长度值可取上文所述六种,主要是考虑嵌入式设备处理能力,一般来讲8位机用256点,这样可以保证下标在8bits以内,16位机可以考虑稍大的点数。
类型数据位数通过:
#define BIT 位数
决定,只能取8,16,32
1.对于图形界面演示示例,信号源和绘图程序在WaveAnalysisView.cpp中的void CWaveAnalysisView::OnDraw(CDC* pDC)函数中。
测试数据采用两个正弦波的叠加,即sin(2*3.14159*i/8)+sin(2*3.14159*i/64);
并将其扩大到最大值之间。
绘图的版本频域图像模值没有开根号,表示的是模值的平方。
如果时域直流量过大,则会引起频率的0点值过大,图像变得缩放的太小,如果为了绘图方便,可将其归零。
2.控制台下示例版本是将timefft.h文件独立出来,和一个控制台结合形成的。 其演示的是一个门函数引发的sinc函数。
定点算法原理和精度分析
由于定点运算位数为定值,因此不能像浮点一样进行各种乘法和加法处理,因为这两种运算会产生溢出。同时,对于旋转因子的小数模式也需要修整为可以运算的整数模式,这便是FFT定点的两个需要处理的难点。
本算法采用了运算量较大但效果较好的动态溢出检查。也就是说在计算每一级之前,对所有值进行遍历,如果发现有超过溢出允许输入值的变量,则对整个数组进行右移缩小处理。这样就保证了数值较小的输入不会被过度缩放,而数值很大的数据也能够通过FFT运算。由于对所有的数组值进行缩放,因此不会造成失真,但是幅度大小会有变化。因此可以理解为定点算法只求解相对值。在使用过程中,尽量不要叠加过多直流分量,因为这将导致频率0点处值很大,其它的值被缩放殆尽。影响精度。
下面详细介绍缩放方法:
一个乘加运算可以理解为:
X1+X2*WN
我们需要找出在什么情况下输出会最大,并在设定输出在此数据类型允许值下时,输入最大取多少。
由于X1,X2,WN都为复数,且abs(X2)=abs(X2*WN),因此最坏的情况是,在模一定的情况下,原来X2的实部和虚部的数值经过旋转因子,被全部转移到了实部,而虚部为零。或者反之。
而由于当X2的实部和虚部相等时,它的模值最大,假设全部转移到实部,在加上X1的实部的值就构成了最大的输出。
我们另这个最大的输出等于类型最大允许值(8bit时为0x7f,16bit时为0x7fff),由于我们需要对所有输入做统一的限定值以便程序判断,所以只能有X1.real=X2.realSqrt(X2.real^2+X2.imag^2)+X1.real= Sqrt(2*X2.real^2)+X1.real=(sqrt(2)+1)X1.real
我们将此输出值限定在类型数据允许的最大整数上,则有
(sqrt(2)+1)X1.real其中(sqrt(2)+1)为常数=2.4142
得到X1.real最大不能超过MAX/2.4142,虚部一样。我们定义此值为OVS
蝶形运算每算完一级,就要进行溢出判断,如果发现其中任何一个值大于此值,则需要对所有值缩放处理。
这时,我们面临怎样缩放的问题,如果采用溢出时所有值/2的方法,那么,它的输入端允许的最大信号幅值就是几乎每级计算完后都需要除以2才通过所有计算。此时
输入/2*2.4142=输出
如果只计算一级,那么输入不得超过OVS/1.207=52(以8bit为例),
而每计算一级都要乘以1.207,成为底数大于1等比数列,此时如果计算1024点的FFT,则要保证输出不溢出,则:
输入*1.207^log(2,N)=输入*1.207^10以八位为例max=127(0x7f),则输入不得超过19,可见,信噪比太低,不能应用。
继续观察运算我们发现,最多可能在同一级上溢出两个bit,我们动态的调整缩放的量级则会改善上述情况。也就是在溢出两位时除以4,溢出一位时除以2,经过计算机模拟,发现此时只要满足单极不溢出的条件(如8bit时为52)则总体就不会溢出,但是,付出的代价是两个bit的精度被处理掉了。因此整体的精度为7bit-2bit(8bit-符号位)=5bit,当然误差可能由于多次这样的操作而累积。
程序中,实现时在每次进入新的一级运算之前,将所有值遍历,如果超出2倍的ovs值则除以4,超出一倍除以2,否则不处理。这样,乘加运算在数据类型内部不会造成溢出了。
旋转因子的缩放:由于在计算旋转因子时,不能用小数,因此,将小数扩大到与位数相当,之后与x[l]相乘,结果再以相同比例缩小,
比如32位版本中:首先将其扩大
wn[k].real=(int)(cos(2*3.14159*k/N)*0x7fffffff);
在计算关于旋转因子的乘积项时乘完的结果再缩小。
XkWn.real=(((__int64)x[l].real*wn[wi].real)>>31)-(((__int64)x[l].imag*wn[wi].imag)>>31);
XkWn.imag=(((__int64)x[l].imag*wn[wi].real)>>31)+(((__int64)x[l].real*wn[wi].imag)>>31);
(__int64)是由于VC不支持32位乘法时自动切换到64为导致了数据丢失,因此做强制转换,其它类型不需要。
matlab中fft定点运算,可用于嵌入式计算的定点FFT算法 (转载)相关推荐
- matlab中fitrsvm函数,训练用于一类和二类分类的支持向量机 (SVM) 分类器
使用 fitcsvm 自动优化超参数. 加载 ionosphere 数据集. load ionosphere 通过使用自动超参数优化,找到最小化五折交叉验证损失的超参数.为了实现可再现性,请设置随机种 ...
- Matlab中fft作频谱横纵坐标
关于这个问题,在很早之前就分享过,也通过了解实现了算法,当时看的明白,想的明白,突然要用的时候,又开始疑问,不免有些纠结,与其每次使用的时候都查,浪费时间,还不如,一次搞定. 真心没把哪门没学好的课程 ...
- matlab变量区表示函数,MATLAB中的工作区,变量和函数
本文概述 工作空间 工作区包含我们在MATLAB中工作时创建的所有变量. 每当我们为变量分配值时, 它都会自动在工作空间中获取空间. 关闭环境后, 工作空间变量将消失, 因此请将这些变量保存在文件中以 ...
- cumsum在matlab中,matlab中cumsum函数和sum函数详解
调用格式及说明 matlab中cumsum函数通常用于计算一个数组各行的累加值.在matlab的命令窗口中输入doc cumsum或者help cumsum即可获得该函数的帮助信息. 调用格式及说明 ...
- MATLAB中颜色模型介绍级各模型之间转换(RGB、HSV、NTSC、YCbCr、HSI)
1.颜色模型定义 2.各颜色模型简介 3.颜色模型的转换 一.颜色模型定义 颜色模型:某个三维颜色空间中的一个可见光子集,它包含某个颜色域的所有颜色.例如,RGB颜色模型就是三维直角坐标颜色系统的一个 ...
- matlab中simple函数怎么用,matlab里simple函数
值 realmin:系统所能表示的最小数值 nargin: 函数的输出引数个数 ---MATLAB 中基本绘图函数有: plot: x 轴和 y 轴均为线性刻度 数刻度 semilogx: x 轴为对 ...
- matlab中sym看不到值和属性,matlab 用sym定义了x,但是输入函数却显示“未定义函数或变量 'x'”?...
答:亲,是syms x,或者是sym('x')来定义x是符号变量 答:matlab2018a中出现未定义函数或"ploy2sym",怎么改要分情况 情况一:符号变量 必须要定义,定 ...
- int在matlab中的作用,int函数表达的是什么意思
int函数表达的是什么意思 int函数相信不少人都没听说过,更别说会知道它表达的是什么意思.为此百分网小编为大家带来int函数表达的意思. int函数表达的意思 C/C++编程语言中,int表示整型变 ...
- matlab 求虚数相位角,在matlab中怎么计算其相位
本文收集整理关于在matlab中怎么计算其相位的相关议题,使用内容导航快速到达. 内容导航: Q1:相位超前补偿器在matlab中是什么模块 首先介绍一下函数,angle()是求相位角,angle() ...
- zeros(0 5)函数matlab,matlab中zeros函数用法
matlab中zeros函数是用于返回一个double类零矩阵,其用法是:1.在命令行窗口中输入"B=zeros(5)",按回车键可生成一个"5*5"的零矩阵: ...
最新文章
- JQ 全选后获取选中的值_为什么在PBI中还需要切片器之三:Excel切片器之度量值切换...
- 16.04linux 安装微信,Ubuntu 16.04安装微信的过程记录
- ADI官方源码快速搭建demo工程验证设计的正确性
- 前端学习(563):干掉block重叠margin重叠
- 【C++深度剖析教程11】C++学习之编写代码实现复数类
- office另存为pdf的加载项_pdf怎么转换成word?打工人必备的丛林法则
- 用了这么多年的PCA可视化竟然是错的!!!
- uboot主循环main_loop
- 让input支持 ctrl v上传粘贴图片? 让input支持QQ截图或剪切板中的图像数据(Java实现保存)...
- 苹果发文谈iPhone SE的核心竞争力,网友:难道不是便宜吗?
- 软件测试人员必备的60个测试工具清单,果断收藏了!
- [fsevents@^1.2.2] optional install error: Package require os(darwin) not compatible with your platfo
- 小米系列手机(包括红米,黑鲨)开启调试模式
- 从NLP任务中文本向量的降维问题,引出LSH(Locality Sensitive Hash 局部敏感哈希)算法及其思想的讨论...
- Windows10当中的混合现实门户怎么使用 超详细讲解 win10混合现实门户怎么用?
- SpringBootTest遇到的问题----Field userMapper in xxx.service.UserService required a bean of type
- 【Grasshopper基础8】电池的序列化与反序列化 Serilization of Grasshopper Component
- JavaOne美国之行–Session篇
- 四、Amlogic A311D 音频回采信号LOOPBACK指南
- 写一个块设备驱动5,6