matlab 倒位序fft程序,[转载]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。
matlab 倒位序fft程序,[转载]MATLAB的一个FFT程序相关推荐
- 用java实现一个计算器程序_1.2第一个java程序——hello world
第一个java程序--hello world 实现一个java程序,主要有三个步骤:1.编写源代码,2.编译源代码,3.运行.java的源代码必须先编译,然后才能由JVM解析执行.所以我们程序员第一步 ...
- PHP 程序员如何设计一个爬虫程序
A8U几年前接过一个项目,类似一个 PHP 爬虫程序,做一个微博舆情分析系统:要爬取新浪微博,用特定关键词搜索中的页面的微博内容. 那是我第一次接触网络爬虫,根本没有思路,也不了解什么Scrap ...
- 小程序多少钱?一个小程序多少钱?
对于需要开发小程序的企业来说,价格是最重要的.现在各大搜索引擎都有很多公司团队推荐小程序开发.那么小程序多少钱?一个小程序多少钱?让我们简单分析两种开发类型. 小程序多少钱?一个小程序多少钱? 一.中 ...
- 怎么创建小程序?如何创建一个小程序呢?
怎么创建小程序?如何创建一个小程序呢?关于怎么创建小程序,是对于一些初次接触小程序或者是想拥有小程序的企业商家来说首先想要搞清楚的问题,其实创建小程序的方法很简单,下面就给大家简单介绍一下怎么创建小程 ...
- 小程序教程1:初识小程序,快速搭建一个小程序项目
前言: 打算整理一个小程序系列的资料.如何快速搭建一个小程序项目,以及一个小程序项目应该有哪些内容. 相关资料: 微信公众平台微信公众平台,给个人.企业和组织提供业务服务与用户管理能力的全新服务平台. ...
- matlab作动态函数曲线图,[转载]Matlab中使用Plot函数动态画图方法总结
本帖最后由 sonictl 于 2012-12-31 12:18 编辑 请删除我 清楚超靠靠靠 没办法,一会儿限制这不能发表,那不能发表的.... [转载]Matlab中使用Plot函数动态画图方法总 ...
- matlab 图像处理 新浪 应变,[转载]Matlab图像处理小结
经常做做图像处理的东西,时间长了,有些函数就忘了,看到网上有人总结,收藏了. 1. 图像和图像数据 缺省情况下,MATLAB将图像中的数据存储为双精度类型(double),64位浮点数,所需存储量很大 ...
- matlab制作动画并保存,[转载]Matlab中动画的实现、制作和保存(转载Matlab技术论坛)...
Matlab中动画的实现.制作和保存(转载Matlab技术论坛) 本帖由Matlab技术论坛原创,原帖参见 http://www.matlabsky.com/thread-592-1-1.html 光 ...
- matlab surf颜色,科学网[转载]matlab中控制mesh或者surf的颜色
科学网[转载]matlab中控制mesh或者surf的颜色 曲面的顏色越发柔順滑腻,del2 (Z) 根據數值微分來產生矩陣 Z 的 discrete Laplacian, Y,surf 及 mesh ...
- matlab中的timer模块,[转载]Matlab中Timer的使用
Matlab中Timer的使用 鉴于Matlab中缺乏多线程机制,使用Timer无疑是一个很重要的工具,Matlab中Timer是一个Java对象. (1)Timer 的定义 t=timer(); 设 ...
最新文章
- jupyter扩展插件Nbextensions的安装、使用
- java 线程分配_Java多线程原子引用分配
- linux程序崩溃时调用链,Linux 获取并分析程序崩溃时的调用堆栈
- NB-IoT在智慧城市中的应用
- 边工作边刷题:70天一遍leetcode: day 6
- 人人学IoT 助学思维导图
- 快速排序-超级详细代码注释!
- cesium获取点击内容信息_Cesium获取鼠标点击位置(PickPosition)
- mac 安装 Xdebug
- 制作自启动的MSDN简体中文版和英文版win98ISO镜像
- win10主题美化(单)
- 微信服务通知消息找回_企业微信消息不提醒怎么办?怎么打开企业微信消息通知?...
- java 场景面试题_Java面试场景整理收录
- 重标极差分析法matlab,重标极差(RS)分析法估计Hurst指数的有效性检验.doc
- win10计算机管理中没有本地用户和组怎么办
- python中reduce函数的作用_python中reduce()函数的使用方法示例
- mysql 1593_Linux中MySQL主从复制中出现1593错误码的低级错误
- 【Jenkins】win 10 / win 11:Jenkins 的下载、安装、部署(Jenkins 2.134 基于 Java 8)
- 自定义Modal组件
- Another Day 超好听的BGM