logo资料库

FFT算法 matlab 实现.doc

发布时间：2022-06-08 发布人：admin 分类：说明书 资料大小：0.29M 资料格式：doc 举报 版权申诉
第1页 / 共9页
第2页 / 共9页
第3页 / 共9页
第4页 / 共9页
第5页 / 共9页
第6页 / 共9页
第7页 / 共9页
第8页 / 共9页
资料共9页，剩余部分请下载后查看
    MATLAB的一个FFT程序
    MATLAB 的一个 FFT 程序 (2011-08-18 17:11:45) 转 载 ▼ 标签： 分类： DSP 数字信号处理 教育 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 程序有： = 512; PointNum 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; + T1 - T2; IMX(lenNum + 2^i + 1) = IMX(lenNum + 1)
    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 位进 行放大处理，数据放大可能会导致溢出，从而使计算错误，而溢出的极限情况是这样：假设
    分享到：
    赞 收藏

    相关推荐

    资料库

    课程资源

    共收录17145份资料，累计13个分类，关注成员有19位，主要包括：PHP,网络管理,网页制作,Java,.Net,数据库,3G/移动开发,C/C++,游戏开发,嵌入式,讲义,软件测试,专业指导

    热门标签

    PHP 网络管理 网页制作 Java .Net 数据库 3G/移动开发 C/C++ 游戏开发 嵌入式 讲义 软件测试 专业指导

    最新资料