MATLAB的一个FFT程序
FFT信号流图:
程序实现是这样:
程序流程如下图:
首先进行位逆转,其实很简单,就是把二进制的位逆转过来:
Matlab的位逆转程序:
function a=bitreverse(Nbit, num)
%Nbit = 4;
%num = 8;
a = 0;
b = bitshift(1,Nbit-1);
for i = 1:Nbit;
if((bitand(num,1)) == 1)
a = bitor(a,b);
end
num = bitshift(num,-1);
b = bitshift(b,-1);
end;
说明:Nbit是逆转位是几位,num是逆转的数即变量。
三个循环,第一个循环是进行N阶的FFT运算
第二个循环其实就是,每一阶FFT的时候,有多少组DFT对象,拿8点来说,第一阶的时候,有4组DFT对象,到了第二阶,就有2组,到了第三,就是最后一阶,只有一组。
第三个循环,其实是在每一组DFT里边,执行多少次蝶形运算!8点DIT FFT来说,第一阶每组有一个蝶形,第二阶每组有2个,第三阶每组有4个蝶形。所以很容易得到三者的关系
i , j, k 三者,反别表示三层循环,然后得出循环次数的关系
stages = log2(PointNum)
i 从 0到stages – 1 !
j 从 0 到 其实就是
k 从0 到
旋转因子W的选择:
因为根据8点DIT-FFT图,从第一阶到最后一阶,可以总结出一个规律:
都是 N是每组蝶形数据个个数,比如第一阶每组有2个元素,N就是2,第二阶每组4个元素,N就是4等。然后x往往都是从0开始到N/2 – 1;
根据旋转因子的性质,其实可以有每阶段每组都是:
蝶形运算设计:
根据信号流图,得出以下算式:
完成了蝶形运算!
全部的matlab程序有:
PointNum = 512;
PointBitNum = 9;
fs = 1024*2;
t = 0:1:PointNum - 1;
%for u = 1:1:PointNum;
sampletab = cos(2*pi*543*t/fs) + cos(2*pi*100*t/fs) + 0.2 + cos(2*pi*857*t/fs) + cos(2*pi*222*t/fs);
%end
zeros(1,PointNum);
sampletab1 = sampletab;
index = 0;
for i = 1:PointNum
k = i - 1
index = bitreverse(PointBitNum,i - 1)
sampletab(i) = sampletab1(index + 1);
end
%sampletab1
%sampletab
REX = sampletab;
IMX = zeros(1,PointNum);
i = 0; ?T Loop for Log2N stages
j = 0; ?T loop for leach sub-DFT
k = 0; ?T Loop for each butterfly
stages = log2(PointNum);
for i = 0 : stages - 1
lenNum = 0;
for j = 0 : 2^(stages - (i + 1)) - 1
for k = 0 : 2^i - 1
R1 = REX(lenNum + 2^i + 1) * cos(2*pi*k*2^(stages - (i + 1))/PointNum);
R2 = IMX(lenNum + 2^i + 1) * sin(2*pi*k*2^(stages - (i + 1))/PointNum);
T1 = REX(lenNum + 2^i + 1) * sin(2*pi*k*2^(stages - (i + 1))/PointNum);
T2 = IMX(lenNum + 2^i + 1) * cos(2*pi*k*2^(stages - (i + 1))/PointNum);
REX(lenNum + 2^i + 1) = REX(lenNum + 1) - R1 - R2;
IMX(lenNum + 2^i + 1) = IMX(lenNum + 1) + T1 - T2;
REX(lenNum + 1) = REX(lenNum + 1) + R1 + R2;
IMX(lenNum + 1) = IMX(lenNum + 1) - T1 + T2 ;
lenNum = lenNum + 1;
endNum = lenNum + 2^i;
end
lenNum = endNum;
end
end
subplot(3,1,1);
fft(sampletab1, PointNum);
x1 = abs(fft(sampletab1, PointNum));
plot([0 : PointNum/2 - 1], x1(1:PointNum/2));
grid on
subplot(3,1,2);
% [REX IMX]
am = sqrt(abs(REX.*REX) + abs(IMX.*IMX));
plot(0:1:PointNum/2 - 1, am(1:PointNum/2));
grid on
subplot(3,1,3);
plot(t, sampletab);
grid on
我还做了与MATLAB原来带有的FFT做比较:
画出的图如下:
第一个是MATLAB自带的FFT函数频谱图
第二个是我自己设计的FFT频谱图
第三个是信号的时域波形
思想已经有了,我以前也改过人家的FFT的C程序但是不是很理解,打算有机会用C语言实现定点FFT,因为在嵌入式上多数用定点FFT,相应的C++版本应该也会写。
下面是网上的一些设计FFT的资料:
N点基-2 FFT算法的实现方法
从图4我们可以总结出对于点数为N=2^L的DFT快速计算方法的流程:
1.对于输入数据序列进行倒位序变换。
该变换的目的是使输出能够得到X(0)~X(N-1)的顺序序列,同样以8点DFT为例,该变换将顺序输入序列x(0)~x(7)变为如图4的x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)序列。其实现方法是:假设顺序输入序列一次村在A(0)~A(N-1)的数组元素中,首先我们将数组下标进行二进制化(例:对于点数为8的序列只需要LOG2(8) = 3位二进制序列表示,序号6就表示为110)。二进制化以后就是将二进制序列进行倒位,倒位的过程就是将原序列从右到左书写一次构成新的序列,例如序号为6的二进制表示为110,倒位后变为了011,即使十进制的3。第三步就是将倒位前和倒位后的序号对应的数据互换。依然以序号6为例,其互换过程如下:
temp = A(6); A(6) = A(3); A(3) = A(6);
实际上考虑到执行效率,如果对于每一次输入的数据都需要这个处理过程是非常浪费时间的。我们可以采用指向指针的指针来实现该过程,或者是采用指针数组来实现该过程。
2.蝶形运算的循环结构。
从图4中我们可以看到对于点数为N = 2^L的fft运算,可以分解为L阶蝶形图级联,每一阶蝶形图内又分为M个蝶形组,每个蝶形组内包含K个蝶形。根据这一点我们就可以构造三阶循环来实现蝶形运算。编程过程需要注意旋转因子与蝶形阶数和蝶形分组内的蝶形个数存在关联。
3.浮点到定点转换需要注意的关键问题
上边的分析都是基于浮点运算来得到的结论,事实上大多数嵌入式系统对浮点运算支持甚微,因此在嵌入式系统中进行离散傅里叶变换一般都应该采用定点方式。对于简单的DFT运算从浮点到定点显得非常容易。根据式(1),假设输入x(n)是经过AD采样的数字序列,AD位数为12位,则输入信号范围为0~4096。为了进行定点运算我们将旋转因子实部虚部同时扩大2^12倍,取整数部分代表旋转因子。之后,我们可以按照(1)式计算,得到的结果与原结果成比例关系,新的结果比原结果的2^12倍。但是,对于使用蝶形运算的fft我们不能采用这种简单的放大旋转因子转为整数计算的方式。因为fft是一个非对称迭代过程,假设我们对旋转因子进行了放大,根据蝶形流图我们可以发现其最终的结果是,不同的输入被放大了不同的倍数,对于第一个输入x(0)永远也不会放大。举一个更加形象的例子,还是以图4为例。从图中可以看出右侧的X(0)可以直接用下式表示:
从上式我们可以看到不同输入项所乘的旋转因子个数(注意这里是个数,就算是wn^0,也被考虑进去了,因为在没有放大时wn^0等于1,放大后所有旋转因子指数模均不为1,因此需要考虑)。这就导致输入不平衡,运算结果不正确。经查阅相关资料,比较妥善的做法是,首先对所有旋转因子都放大2^Q倍,Q必须要大于等于L,以保证不同旋转因子的差异化。旋转因子放大,为了保证其模为1,在每一次蝶形运算的乘积运算中我们需要将结果右移Q位来抵消这个放大,从而得到正确的结果。之所以采用放大倍数必须是2的整数次幂的原因也在于此,我们之后可以通过简单的右移位运算将之前的放大抵消,而右移位又代替了除法运算,大大节省了时间。
4.计算过程中的溢出问题
最后需要注意的一个问题就是计算过程中的溢出问题。在实际应用中,AD虽然有12位的位宽,但是采样得到的信号可能较小,例如可能在0~8之间波动,也就是说实际可能只有3位的情况。这种情况下为了在计算过程中不丢失信息,一般都需要先将输入数据左移P位进行放大处理,数据放大可能会导致溢出,从而使计算错误,而溢出的极限情况是这样:假设我们数据位宽为D位(不包括符号位),AD采样位数B位,数字放大倍数P位,旋转因此放大倍数Q位,FFT级联运算带来的最大累加倍数L位。我们得到:
假设AD位宽12,数据位宽32,符号位1位,因此有效位宽31位,采样点数N,那么我们可以得到log2(N)+P+Q<=19,假设点数128,又Q>=L可以得到放大倍数P<=5。
转载于:https://www.cnblogs.com/nickchan/archive/2011/08/18/3104549.html
MATLAB的一个FFT程序相关推荐
- matlab程序 如何使用,如何使用MATLAB创建一个最简单的程序
<如何使用MATLAB创建一个最简单的程序>由会员分享,可在线阅读,更多相关<如何使用MATLAB创建一个最简单的程序(4页珍藏版)>请在人人文库网上搜索. 1.如何使用MAT ...
- 用MATLAB写一个自动生成福利彩票双色球号码的程序
用MATLAB写一个自动生成福利彩票双色球号码的程序 规则 红色球:1-33号任选6个 蓝色球:1-16号任选1个 red = randi([1,33],1,6); disp('红色球'); fpri ...
- 用matlab设计一个简单的抽奖程序
问题描述 国庆节快要到了,实验室要求设计一个抽奖程序.先将所有实验室成员的姓名输入到一个excel文件中,然后运行程序,读取excel文件中的数据,从中随机抽取一个人作为中奖者.用matlab就可以解 ...
- 用matlab写一个利用人工地震波定位掩埋物的程序
利用人工地震波定位掩埋物是地质勘探中常用的方法之一.下面是一个使用MATLAB编写的基本程序框架,可以用于实现这一目的. 程序主要分为以下几个步骤: 1.读取地震数据文件 使用MATLAB中的load ...
- matlab编一个福利彩票电脑选号的程序,【Matlab编程】Matlab让电脑失而复得
在学校常常有同学电脑失窃,大抵都是粗细大意.据说iPhone手机失窃后能够获取小偷的照片,从而将照片找到.如今用matlab写一个程序使得当小偷使用电脑上网时,电脑自己主动将电脑前面的人的照片发到你指 ...
- matlab相关性分析频谱_基于Matlab的相关频谱分析程序教程
基于Matlab的相关频谱分析程序教程 Matlab 信号处理工具箱 谱估计专题 频谱分析 Spectral estimation(谱估计)的目标是基于一个有限的数据集合描述一个信号的功率(在频率上的 ...
- matlab有意思程序,matlab有意思的小程序
10个C++趣味小程序,很有意思的.VIP专享文档 VIP专享文档是百度文库认... 现在很多人使用微信的时间已经非常长了,他们注册的微信号往上可能已经是5年前的事情了,正是由于不少使用者在这个过程当 ...
- 车牌识别与计算机编程,基于MATLAB的车牌识别程序详解.ppt
基于MATLAB的车牌识别程序详解 自定义一个字符函数,用来从车牌区域中提取出7个字符,其中利用切割函数来进行切割. 程序:function [word,result]=getword(d) word ...
- matlab控制算法C语言,PID算法Matlab仿真程序和C程序
<PID算法Matlab仿真程序和C程序>由会员分享,可在线阅读,更多相关<PID算法Matlab仿真程序和C程序(6页珍藏版)>请在人人文库网上搜索. 1.增量式PID控制算 ...
- matlab 识别调试,有关matlab的人脸识别程序,但调试是不成功
有关matlab的人脸识别程序,但调试是不成功,求高手帮忙指点修改.在此先谢了 1.色彩空间转换 function[r,g]=rgb_RGB(Ori_Face) R=Ori_Face(:,:,1); ...
最新文章
- golang odbc mysql_golang使用odbc链接hive
- HTML5学习笔记(二十七):Ajax
- SQL Server 2014 Win7 Win10 安装详解 SQL Server 2017 2019 Linux及SQL TSQL ETL实用案例
- Activiti绩效对决
- method=post 怎么让查看源代码看不到_网上文档无法复制怎么办?试试这几个方法!...
- 2009 CCTV体坛风云人物颁奖盛典,精彩语录
- 【Vue.js源码解析 三】-- 模板编译和组件化
- 自制一个简单的操作系统
- ibm x60 学习linux,IBM X60 T60系列安装系统时SATA设置问题
- Linux/软件 - 资源[国外站点]
- 一些有趣但少有人知的 Python 特性
- 华为Mate40/华为Mate40Pro忘记密码怎么解锁激活手机设备已锁定恢复出厂无法解锁账户ID屏幕锁解除刷机方法教程
- ByteBuffer Converting CharBuffer
- Linux内存管理二(页表)
- python中justify的意思_Python3 tkinter基础 Label justify 多行字符串左对齐
- 百度飞桨第一天学习笔记
- 沉浸式过山车:穿梭时空,自在畅游巨蚁数字虚实世界
- 蓝牙5.0对比4.2的主要优势
- 肖特基二极管型号大全之ASEMI肖特基常见型号
- 陌陌、BIGO、比亚迪、好未来、同程艺龙、去哪儿、联通西安研究院、华为、海康威视java开发面试经验
热门文章
- mysql with roll up_GROUP BY...WITH ROLL UP 分组统计后的再合计
- bootdo框架首页解析
- python网页爬虫菜鸟教程_Python爬虫实践(7)-抓取菜鸟教程python学习路线-工具-站长头条...
- pyGurobi使用手册
- ac1900 linksys 恢复_把变砖的Linksys-AC1900路由器救活
- 百度之星冠军分享:AI图像赛事入门
- 汇编语言debug的使用
- linux simhei 字体下载,Linux CentOS 7 安装 字体库文件(simsun.ttf、simheittf.ttf)
- 定时器 cron 表达式
- 2021~ 你好,加油 (ง •_•)ง