离散傅里叶变换及matlab实现(按时间抽选(DIT)的基-2 FFT算法(库利-图基算法))
转,傅里叶变换,很好的解释 很好的文章,可惜水平太差,还没有完全理解。
快速傅里叶的matlab实现
按时间抽选(DIT)的基-2 FFT算法(库利-图基算法)
傅里叶要用到的nn个复数,不是随机找的,而是——把单位圆(圆心为原点、1为半径的圆)nn等分,取这nn个点(或点表示的向量)所表示的虚数,即分别以这nn个点的横坐标为实部、纵坐标为虚部,所构成的虚数。
64点傅里叶变换的matlab代码
function FFT(x)
%******************************************************
%*********按时间抽选的基2-FFT算法*******
%输入参数:x-离散时间信号
%输出参数:X-序列x的N点DFT(N是序列长度,必须是2的整数次幂)
%*******************************************************N = 64; %64点FFT N应该为一个全局变量M=log2(N);if length(x) < N x = [x; zeros(1,N-length(x))']; %若x的长度不是2的整数次幂,则补零至NendBRTable = bin2dec(fliplr(dec2bin([1:N]-1,M))) + 1; %求 1:N 序列序号的倒位序:1.将系数取二进制并对调,然后取十进制 %2.matlab的下标从1开始,需要加1X = x(BRTable); %调整x输入顺序后的序列,并作为X的初始化 WN = exp(-j*2*pi/N); %旋转因子for L = 1:M;B = 2^(L-1); %第L级中,每个蝶形的两个输入数据相聚B个点,共有B个旋转因子for J = 0:B-1; %第L级中不同的旋转因子p = J*2^(M-L); %旋转因子的指数WNp = WN^p;for k = J+1 : 2^L : N ; t = X(k+B) * WNp; %产生错误:未定义与'cell'类型的输入参数相对应的运算符‘*’-x=[]非x={}X(k+B) = X(k)-t;X(k) = X(k)+t; endendend
end
离散化处理,求64点FFT(以下代码已更新,之前有几个符号写错了,天知道我是怎么工作的!!!以头抢地尔)
function Result = FFT(x, N)
x = double(x);
f = 1;
if f > 0%N = 64; %64点FFT N应该为一个全局变量M = log2(N);if length(x) < N ;x = [x zeros(1,N-length(x))]; %若x的长度不是2的整数次幂,则补零至NendBRTable = bin2dec(fliplr(dec2bin([1:N]-1,M))) + 1; %求 1:N 序列序号的倒位序:1.将系数取二进制并对调,然后取十进制 %2.matlab的下标从1开始,需要加1Xreal = x(BRTable); %调整x输入顺序后的序列,并作为X的初始化 Ximag = zeros(1,N);X = x(BRTable);WN = exp(-j*2*pi/N);for L = 1:M;B = 2^(L-1); %第L级中,每个蝶形的两个输入数据相聚B个点,共有B个旋转因子for J = 0:B-1; %第L级中不同的旋转因子p = J*2^(M-L); %旋转因子的指数WNp = WN^p;%p = p*j;for k = J+1: 2^L : N; %计算出来的包括实部和虚部 此种方法已验证正确f=0;if f>0t = X(k+B) * WNp; %产生错误:未定义与'cell'类型的输入参数相对应的运算符‘*’-x=[]非x={} X(k+B) = X(k)-t; %将公式写成指数的形式,按虚实区分X(k) = X(k)+t;else %计算出来的实部 虚部分开,正确Wncosp = cos(2*pi/N*p);Wnsinp = sin(2*pi/N*p);Treal = Xreal(k+B)*Wncosp+Ximag(k+B)*Wnsinp;Timag = Ximag(k+B)*Wncosp-Xreal(k+B)*Wnsinp;Xreal(k+B) = Xreal(k) -Treal;Ximag(k+B) = Ximag(k) -Timag;Xreal(k) = Xreal(k)+Treal;Ximag(k) = Ximag(k)+Timag;%以下C语言中应该可以
% Xreal(k+B) = Xreal(k)-Xreal(k+B)*cos(2*pi/N*p)-Ximag(k+B)*sin(2*pi/N*p);
% Ximag(k+B) = Ximag(k)+Xreal(k+B)*sin(2*pi/N*p)-Ximag(k+B)*cos(2*pi/N*p);
% Xreal(k) = Xreal(k)+Xreal(k+B)*cos(2*pi/N*p)+Ximag(k+B)*sin(2*pi/N*p);
% Ximag(k) = Ximag(k)-Xreal(k+B)*sin(2*pi/N*p)+Ximag(k+B)*cos(2*pi/N*p);%将l写成了1,以下正确
% a1 = real(X(k+B));
% a2 = real(exp(-j*2*pi/N*p));
% b1 = imag(X(k+B));
% b2 = imag(exp(-j*2*pi/N*p));
% temp =a1*a2-b1*b2+j*(a1*b2+b1*a2);
% X(k+B) = X(k)-temp;
% X(k) = X(k)+temp;
endendendend%}
elseL = ceil(log2(length(x)));N = 2^L;if length(x) < Nx = [x zeros(1,N-length(x))];end BRTable = bin2dec(fliplr(dec2bin([1:N]-1,L))) + 1;% x_rep = x(BRTable);x_rep = seq_replace(x);for l = 1:Ld = 2^(l-1);for n = 1:dfor k = n:2*d:Nk_dual = k+d;r = (k-1)*2^(L-1);a1 = real(x_rep(k_dual));a2 = real(exp(-j*2*pi/N*r));b1 = imag(x_rep(k_dual));b2 = imag(exp(-j*2*pi/N*r));temp = a1*a2-b1*b2+j*(a1*b2+b1*a2);x_rep(k_dual) = x_rep(k) -temp;x_rep(k) = x_rep(k)+temp;endendendend%求数据的模for i = 1:1:N;Rereal(i) = Xreal(i)*Xreal(i);Reimag(i) = Ximag(i)*Ximag(i);Result(i) = sqrt(Rereal(i)+Reimag(i));%fprintf('%d ',Result(i));end
end
输入错一个:,写成了;,找了一小时的错误,唉
基于64点输入如下:
L = log2(N);if length(x) < Nx = [x zeros(1,N-length(x))];
end BRTable = bin2dec(fliplr(dec2bin([1:N]-1,L))) + 1;
x_rep = x(BRTable);for l = 1:Ld = 2^(l-1);for n = 1:dfor k = n:2*d:Nk_dual = k+d;r = (k-1)*2^(L-l); %这里是l,不是1a1 = real(x_rep(k_dual));a2 = real(exp(-j*2*pi/N*r));b1 = imag(x_rep(k_dual));b2 = imag(exp(-j*2*pi/N*r));temp = a1*a2-b1*b2+j*(a1*b2+b1*a2);x_rep(k_dual) = x_rep(k) -temp;x_rep(k) = x_rep(k)+temp;endend
end
验证
x=[1 2 3 4 5 6 7 8];
调用 FFT(x,8);
result =1 至 7 列36.0000 + 0.0000i -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i -4.0000 + 0.0000i -4.0000 - 1.6569i -4.0000 - 4.0000i8 列-4.0000 - 9.6569i
基于64点的FFT求128点FFT
%% 128点FFT(基于以上64点)
% x:输入序列
% N:128
% X:输出序列
function X = FFT128(x,N)if length(x) < Nx = [x,zeros(1,N-length(x))];endy1 = FFT(x(1:2:128),64); %计算偶数组的64点FFTy2 = FFT(X(2:2:128),64); %计算奇数组的64点FFTwn = exp(-j*2*pi/128);X(1:64) = y1+y2.*(wn.^[0:63]);X(65:128) = y1-y2.*(wn.^[0:63]);end
1.
2,
3,
离散傅里叶变换及matlab实现(按时间抽选(DIT)的基-2 FFT算法(库利-图基算法))相关推荐
- 按时间抽选(DIT)的基-2 FFT算法(库利-图基算法)C++程序
基-2 FFT算法的C++程序,按时间抽选.输入倒位序.输出自然顺序,N=2LN=2^LN=2L #include <complex>int fft(complex<double&g ...
- 2.MATLAB利用“基2时间抽选法”实现FFT
通过基2时间抽选法的原理,编程实现基2(DIT)FFT 文章目录 题目重述 问题分析以及求解思路 程序代码 题目重述 问题分析以及求解思路 待完善(请耐心等待) 程序代码 %%数据倒位序 N=32; ...
- c语言编程实现基2-fft,时间抽选基2FFT及IFFT算法C语言实现
/*时间抽选基2FFT及IFFT算法C语言实现*/ /*Author :Junyi Sun*/ /*Copyright 2004-2005*/ /*Mail:ccnusjy@yahoo.com.cn* ...
- c代码实现 ifft运算_月光软件站 - 编程文档 - 其他语言 - 时间抽选基2FFT及IFFT算法C语言实现...
正在学数字信号处理,感觉上学期信号与系统学得不扎实,因为当时只是死记公式,这学期数信老师提倡动手实践,觉得自己在编程中对公式理解得更加深刻了. 以下是我写的FFT,欢迎指教. /*时间抽选基2FFT及 ...
- matlab基2时间抽选法,按时间抽取的基2FFT算法分析及MATLAB实现
电子技术研发ElectronicsR&D 电一子一技一术- 按时问抽取的基2FFT算法分析及MATLAB实现 张登奇李宏民李丹 (湖南理工学院信息与通信工程学院) 摘要:DFT是一种应用广泛的 ...
- 用MATLAB计算序列的离散傅里叶变换
用MATLAB计算序列的离散傅里叶变换 MATLAB提供了用快速算法计算离散傅里叶变换的函数fft,其调用格式为: Xk = fft(xn, N) 其中,调用参数xn为时域序列向量,N为离散傅里叶变换 ...
- 图像处理-离散傅里叶变换-数字图像处理第三版第四章内容
图像傅里叶变换方法有很多,可以通过空间光调制器输入图像后在通过平行光照明经过傅里叶变换透镜进行傅里叶变换,另一个方法就是利用计算机进行傅里叶变换,其中傅里叶变换有两种算法一种是DFT还有一种是FFT( ...
- 【从FT到DFT和FFT】(三)从离散傅里叶变换到快速傅里叶变换
文章目录 推荐阅读 前言 从离散傅里叶变换到快速傅里叶变换 单位根 对DFT进行分治得到FFT 计算前半截 计算后半截 快速傅里叶逆变换(IFFT) 推荐阅读 前置阅读 [从FT到DFT和FFT](一 ...
- Matlab如何进行利用离散傅里叶变换DFT (快速傅里叶变换FFT)进行频谱分析
文章目录 1. 定义 2. 变换和处理 3. 函数 4. 实例演示 例1:单频正弦信号(整数周期采样) 例2:单频正弦信号(非整数周期采样) 例3:含有直流分量的单频正弦信号 例4:正弦复合信号 例5 ...
最新文章
- [转] ASP.NET MVC3 路由和多数据集的返回
- 转:去掉Flex4生成的SWF加载时的进度条
- jQuery Mobile 图标无法显示
- 【STM32】FreeRTOS 中断配置和临界段
- React开发(239):dva概念4dispatch
- 信息学奥赛C++语言:数字反转
- 给网站logo添加css帅气亮光扫过特效 附教程
- shell-script(command groups)
- 360解压电脑版安装包_迅捷pdf转换器电脑版安装包下载-迅迅捷pdf转换器安装包免费下载...
- UVA11152 Safe Salutations【计算几何】
- 投屏后能在电脑操作手机吗 手机投屏电脑操作手机软件
- sap销售发货的流程_现金及快速销售流程
- python群发邮件
- 飞秋2012、飞秋2013资源文件
- 从链家网上爬取租房数据并进行数据分析
- 网站采集器-免费任意网页数据采集器
- php制作万年历的步骤_PHP制作万年历
- Socket编程中具体接口的用法
- Linux 命令(217)—— iptables-restore 命令
- 「硬见小百科」几款常用的保护电路