快速傅里叶变换_蝶形运算_按时间抽取基2-fft算法_MATLAB代码
- 格式:doc
- 大小:36.50 KB
- 文档页数:4
离散时间信号的基2快速傅里叶变换FFT (时间抽取)蝶形算法实现一、一维连续信号的傅里叶变换连续函数f(x)满足Dirichlet (狄利克雷)条件下,存在积分变换:正变换:2()()()()j ux F u f x e dx R u jI u π+∞--∞==+⎰ 反变换:2()()j ux f x F u e du π+∞-∞=⎰其中()()cos(2)R u f t ut dt π+∞-∞=⎰,()()sin(2)I u f t ut dt π+∞-∞=-⎰定义幅值、相位和能量如下:幅度:1222()()()F u R u I u ⎡⎤⎡⎤=+⎣⎦⎣⎦ 相位:()arctan(()/())u I u R u ϕ= 能量:22()()(E u R u I u =+)二、一维离散信号的傅里叶变换将连续信号对自变量进行抽样得到离散信号(理想冲击抽样脉冲),利用连续信号的傅里叶变换公式得到离散时间傅里叶变换DTFT ;再利用周期性进行频域抽样,得离散傅里叶变换DFT (详情参考任何一本《数字信号处理》教材)。
DFT 变换如下:正变换:12/0()(),0,1,2,1N j ux Nx F u f x eu N π--===-∑。
反变换:12/01()(),0,1,2,1N j ux Nu f x F u ex N Nπ-===-∑。
DFT 是信号分析与处理中的一种重要变换,因为计算机等数字设备只能存储和处理离散数据(时域、频域)。
因直接计算DFT 的计算量大(与变换区间长度N 的平方成正比,当N 较大时,计算量太大),所以在快速傅里叶变换(简称FFT)出现以前,直接用DFT 算法进行谱分析和信号的实时处理是不切实际的。
直到1965年发现了DFT 的一种快速算法(快速傅里叶变换,即FFT )以后,情况才发生了根本的变化。
FFT 有时间抽取和频率抽取两种,下面介绍时间抽取FFT 。
三、时间抽取的基2快速傅里叶变换FFT令2j NN W eπ-=,则2jkm km NNWeπ-=称为旋转因子,把DFT 正变换改写为:1[][],0,1,1N km N k F m f k W m N -===-∑将函数记作x ,变换后为X ,则为:10[][],0,1,1N kmN k X m x k W m N -===-∑时间抽取的FFT 算法利用了旋转因子的三个性质:周期性、对称性和可约性。
1 概述2 代码3 算例1 概述任何连续测量的时序或信号,都可以表示为不同频率的余弦(或正弦)波信号的无限叠加。
FFT (Fast Fourier Transform )是离散傅立叶变换的快速算法,可以将一个信号变换到频域。
对于包含 个均匀采样点的向量 ,其傅里叶变换定义为式中:,为虚数单位为什么做FFT :(1)有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征(频率,幅值,初相位);(2)FFT 可以将一个信号的频谱提取出来,进行频谱分析,为后续滤波准备;(3)通过对一个系统的输入信号和输出信号进行快速傅里叶变换后,两者进行对比,对系统可以有一个初步认识。
2 代码function [ExtractedSignal]=FFTransform(OriginalSignal,T,Frequency,varargin)% description:% [ExtractedSignal]=FFTransform(OriginalSignal,T,Frequency,Threshold)% 快速傅里叶变换提取信号% input:% OriginalSignal-----观测值序列% T------------------采样间隔% Frequency----------欲提取的信号频率,向量% varargin-----------可选参数Threshold ,频率阈值,默认为1e-6,% 原始信号频率与给定频率差值大于该阈值则予以剔除% output:% ExtractedSignal-----FFT 变换后提取的信号%%p =inputParser;addOptional(p,'Threshold',1e-6);parse(p,varargin{:});Threshold =p.Results.Threshold;12345678910111213141516171819203 算例假设一个随时间t 变化的信号。
按时间抽取的基2FFT算法分析基2FFT算法是一种快速傅里叶变换算法,它通过将傅里叶变换的计算复杂度从O(n^2)降低到O(nlogn),大大提高了傅里叶变换的效率。
基2FFT算法的核心思想是将一个长度为n的序列分成长度为n/2的两个子序列,并分别做傅里叶变换。
然后将两个子序列的傅里叶变换结果合并起来,得到原始序列的傅里叶变换结果。
具体来说,基2FFT算法的步骤如下:1.如果输入序列长度为1,则返回输入序列作为傅里叶变换结果。
2.将输入序列按奇偶位置分为两个子序列。
3.对两个子序列分别递归地应用基2FFT算法,得到它们的傅里叶变换结果。
4.根据蝶形算法,将子序列的傅里叶变换结果合并起来,得到原始序列的傅里叶变换结果。
基2FFT算法通过不断将序列分成两半的方式,将傅里叶变换的计算复杂度从O(n^2)降低到O(nlogn)。
在每一层递归中,需要进行O(n)次计算,而递归的层数为logn,因此总的时间复杂度为O(nlogn)。
基2FFT算法的关键之一是蝶形算法。
蝶形算法是一种合并子序列傅里叶变换结果的方法。
在每一层递归中,对于每个位置k,需要计算一个长度为n的序列上的k点DFT。
根据蝶形算法,可以将这个计算分成两个部分:计算序列中偶数位置上的点DFT和计算序列中奇数位置上的点DFT,并通过一些乘法和加法操作合并起来。
这样做可以大大减少计算量,提高计算效率。
基2FFT算法还可以通过多线程或并行处理来进一步提高效率。
由于基2FFT算法具有递归结构,可以将不同的递归层分配给不同的线程或处理器来并行进行计算,从而加快计算速度。
基2FFT算法在数字信号处理、图像处理、通信系统和科学计算等领域有着广泛的应用。
它的高效性和快速运算速度使得它成为处理大规模数据的重要工具。
综上所述,基2FFT算法通过将傅里叶变换的计算复杂度从O(n^2)降低到O(nlogn),大大提高了傅里叶变换的效率。
它采用递归分治的思想,通过分解和合并操作来实现傅里叶变换的计算。
一、实验目的1.利用MATLAB 编写FFT 快速傅里叶变换。
2.比较编写的myfft 程序运算结果与MATLAB 中的FFT 的有无误差。
二、实验条件PC 机,MATLAB7.0三、实验原理1. FFT (快速傅里叶变换)原理:将一个N 点的计算分解为两个N/2点的计算,每个N/2点的计算再进一步分解为N/4点的计算,以此类推。
根据DFT 的定义式,将信号x[n]根据采样号n 分解为偶采样点和奇采样点。
设偶采样序列为y[n]=x[2n],奇采样序列为z[n]=x[2n+1]。
上式中的k N W -为旋转因子N k j e /2π-。
下式则为y[n]与z[n]的表达式:2.蝶形变换的原理:下图给出了蝶形变换的运算流图,可由两个N/2点的FFT(Y[k]和Z[k]得出N点FFT X[k])。
同理,每个N/2点的FFT可以由两个N/4点的FFT求得。
按这种方法,该过程可延迟后推到2点的FFT。
下图为N=8的分解过程。
图中最右边的为8个时域采样点的8点FFTX[k],由偶编号采样点的4点FFT和奇编号采样点的4点得到。
这4点偶编号又由偶编号的偶采样点的2点FFT和奇编号的偶采样点的2点FFT产生。
相同的4点奇编号也是如此。
依次往左都可以用相同的方法算出,最后由偶编号的奇采样点和奇编号的偶采样点的2点FFT算出。
图中没2点FFT成为蝶形,第一级需要每组一个蝶形的4组,第二级有每组两个蝶形的两组,最后一级需要一组4个蝶形。
四、实验内容1.定义函数disbutterfly ,程序根据FFT 的定义:]2[][][N n x n x n y ++=、n N W N n x n x n z -+-=])2[][(][,将序列x 分解为偶采样点y 和奇采样点z 。
function [y,z]=disbutterfly(x)N=length(x);n=0:N/2-1;w=exp(-2*1i*pi/N).^n;x1=x(n+1);x2=x(n+1+N/2);y=x1+x2;z=(x1-x2).*w;2.定义函数rader ,纠正输出序列的输出顺序。
快速傅里叶变换FFT傅里叶变换是频谱分析、信号分析、线性系统分析的重要工具。
傅里叶变换在物理学、电子类学科、信号处理、声学、结构动力学等领域都有着广泛的应用。
将时域信息写成Fourier 级数:01()(cos sin )2T n n n a f t a n t b n t ωω∞==++∑ 222222022()d ,,2()cos (1,2,),2()sin (1,2,).T T TT TT T n T n T a f t t T Ta f t n t dt n Tb f t n t dt n T πωωω---======⎰⎰⎰其中 0001(1)[(0)(0)].2T T t f t f t ++-在间断点处,式右端级数收敛于引入欧拉公式=cos sin cos ,sin .22i i i i i e i e e e e i φφφφφφφφφ--++-==- 00,2,,1,2,3,,22n n n n n n a c a i b a i b c c n -=-+===令 则()in t T n n f t c e ω+∞=-∞=∑Fourier 级数的复指数形式具有明显的物理意义。
离散傅里叶变换DFT为数值计算n c ,需要进行离散,离散后的形式为()()()012/011N a T jn N k jn k t n a n c k t e dt k t e t T N t πωωω-+--⋅∆=≈∆=⋅∆⋅∆⋅∆∑⎰记()()()12/0N jnk N n X n k t eπω--==⋅∆∑考虑到()X n 本身以N 为周期,于是复数形式Fourier 级数展开式中()()102=102n N X n n N c N X N n n N⎧≤<⎪⎪⎨⎪+-≤<⎪⎩ 变换后的幅值/2n n C A == 真实幅值()22n n A C X n N==快速傅里叶变换FFTDFT 计算量大,为了降低计算量,在DFT 的基础上利用系数的对称性和周期性,将长序DFT 转为短序DFT 。
实验14 快速傅里叶变换(FFT)(完美格式版,本人自己完成,所有语句正确,不排除极个别错误,特别适用于山大,勿用冰点等工具下载,否则下载之后的word 格式会让很多部分格式错误,谢谢)XXXX 学号姓名处XXXX一、实验目的1、加深对双线性变换法设计IIR 数字滤波器基本方法的了解。
2、掌握用双线性变换法设计数字低通、高通、带通、带阻滤波器的方法。
3、了解MA TLAB 有关双线性变换法的子函数。
二、实验内容1、双线性变换法的基本知识2、用双线性变换法设计IIR 数字低通滤波器3、用双线性变换法设计IIR 数字高通滤波器4、用双线性变换法设计IIR 数字带通滤波器三、实验环境MA TLAB7.0四、实验原理1、实验涉及的MATLAB 子函数(1)fft功能:一维快速傅里叶变换(FFT)。
调用格式:)(x fft y =;利用FFT 算法计算矢量x 的离散傅里叶变换,当x 为矩阵时,y 为矩阵x每一列的FFT 。
当x 的长度为2的幂次方时,则fft 函数采用基2的FFT 算法,否则采用稍慢的混合基算法。
),(n x fft y =;采用n 点FFT 。
当x 的长度小于n 时,fft 函数在x 的尾部补零,以构成n点数据;当x 的长度大于n 时,fft 函数会截断序列x 。
当x 为矩阵时,fft 函数按类似的方式处理列长度。
(2)ifft功能:一维快速傅里叶逆变换(IFFT)。
调用格式:)(x ifft y =;用于计算矢量x 的IFFT 。
当x 为矩阵时,计算所得的y 为矩阵x 中每一列的IFFT 。
),(n x ifft y =;采用n 点IFFT 。
当length(x)<n 时,在x 中补零;当length(x)>n 时,将x 截断,使length(x)=n 。
(3)fftshift功能:对fft 的输出进行重新排列,将零频分量移到频谱的中心。
调用格式:)(x fftshift y =;对fft 的输出进行重新排列,将零频分量移到频谱的中心。
7.4实验4:快速傅里叶变换-基2时间抽取FFT 算法matlab 实现7.4.1实验目的1.练习利用matlab6.5中工具箱中的信号处理函数2.熟悉快速傅里叶变换的基本原理3.熟悉基2DIT-FFT 运算的MATLAB 程序并运用7.4.2涉及函数信号处理函数X=fft(x)或者X=fft(x,N):自定义功能函数function y=myfft(xr,n)7.4.3实验原理与方法1 DIT-FFT 算法的基本原理有限长序列x (n )的N 点DFT 定义为:∑-==10 )()(N n n k N W n x k X ,式中j N eW π2-=,其整数次幂简称为旋转因子。
N 符合2的整数幂,N 为2的几次幂,则需要进行几次分解。
碟形运算流图符号如下:2 DIT-FFT 算法的运算规律及编程思想为了编写DIT-FFT 算法的运算程序,首先要分析其运算规律,总结编程思想并绘出程序框图。
由右图可知,DIT-FFT 算法的运算过程很有规律。
2.1 原位计算对M N 2=点的FFT 共进行M 级运算,每级由N /2个蝶形运算组成。
在同一级中,每个蝶的输入数据只对本蝶有用,且输出节点与输入节点在同一水平线上,这就意味着每算完一个蝶后,所得数据可立即存入原输入数据所占用的数组元素(存储单元),这种原位(址)计算的方法可节省大量内存。
2.2 蝶形运算实现FFT 运算的核心是蝶形运算,找出蝶形运算的规律是编程的基础。
for mm=1:m %将DFT 做m 次基2分解,从左到右,对每次分解作DFT 运算 Nmr=2^mm;u=1; %旋转因子u 初始化WN=exp(-j*2*pi/Nmr); %本次分解的基本DFT 因子WN =exp(-i*2*pi/Nmr)for n=1:Nmr/2 %本次跨越间隔内的各次碟形运算for k=n:Nmr:N %本次碟形运算的跨越间隔为Nmr=2^mmkp=k+Nmr/2; %确定碟形运算的对应单元下标(对称性)t=x(kp)*u; %碟形运算的乘积项x(kp)=x(k)-t; %碟形运算的加法项x(k)=x(k)+t;endu=u*WN; %修改旋转因子,多乘一个基本DFT 因子WN2.3 序列倒序为了保证运算输出的X (k )按顺序排列,要求序列x (n )倒序输入,即在运算前要先对输入的序列进行位序颠倒。
按时间抽取的基2FFT算法分析及MATLAB实现基2FFT算法是一种快速傅里叶变换(Fast Fourier Transform,FFT)的算法,在信号处理、图像处理等领域有着广泛的应用。
该算法通过将N个输入值分解成两个长度为N/2的DFT(离散傅里叶变换)来实现快速的计算。
本文将对基2FFT算法进行分析,并给出MATLAB实现。
基2FFT算法的主要思路是将输入序列分解成奇偶两个子序列,然后分别对这两个子序列进行计算。
具体步骤如下:1.将输入序列拆分成奇数位和偶数位两个子序列。
比如序列x[0],x[1],x[2],x[3]可以拆分成x[0],x[2]和x[1],x[3]两个子序列。
2. 对两个子序列分别进行DFT计算。
DFT的定义为:X[k] = Σ(x[n] * exp(-i * 2π * k * n / N)),其中k为频率的索引,N为序列长度。
3.对得到的两个DFT结果分别进行合并。
将奇数位子序列的DFT结果和偶数位子序列的DFT结果合并成一个长度为N的DFT结果。
4.重复以上步骤,直到计算结束。
基2FFT算法的时间复杂度为O(NlogN),远远小于直接计算DFT的时间复杂度O(N^2)。
这是因为基2FFT算法将问题的规模逐步减半,从而实现了快速的计算。
下面是MATLAB中基2FFT算法的简单实现:```matlabfunction X = myFFT(x)N = length(x);if N == 1X=x;%递归结束条件return;endeven = myFFT(x(1:2:N)); % 偶数位子序列的FFT计算odd = myFFT(x(2:2:N)); % 奇数位子序列的FFT计算W = exp(-1i * 2 * pi / N * (0:N/2-1)); % 蝶形因子temp = W .* odd; % 奇数位子序列的DFT结果乘以蝶形因子X = [even + temp, even - temp]; % 合并得到一个长度为N的DFT结果end```上述代码中,函数myFFT为基2FFT算法的MATLAB实现。
function y=MyFFT_TB(x,n)
%MYFFT_TB:My Fast Fourier Transform Time Based
%按时间抽取基2-fft算法
%input:
% x -- 输入的一维样本
% n -- 变换长度,缺省时n=length(x) 当n小于x数据长度时,x数据被截断到第n个数据% 当n大于时,x数据在尾部补0直到x 含n个数据
%output:
% y -- 1*n的向量,快速傅里叶变换结果
%variable define:
% N -- 一维数据x的长度
% xtem -- 临时储存x数据用
% m,M -- 对N进行分解N=2^m*M,M为不能被2整除的整数
% two_m -- 2^m
% adr -- 变址,1*N的向量
% l -- 当前蝶形运算的级数
% W -- 长为N/2的向量,记录W(0,N),W(1,N),...W(N/2-1,N)
% d -- 蝶形运算两点间距离
% t -- 第l级蝶形运算含有的奇偶数组的个数
% mul -- 标量,乘数
% ind1,ind2 -- 标量,下标
% tem -- 标量,用于临时储存
%参考文献:
% /view/fea1e985b9d528ea81c779ee.html
%% 输入参数个数检查
msg=nargchk(1,2,nargin);
error(msg);
%% 输入数据截断或加0
N=length(x);
if nargin==2
if N<n % 加0
xtem=x;
x=zeros(1,n);
x(1:N)=xtem;
N=n;
else % 截断
xtem=x;
x=xtem(1:n);
N=n;
end
end
%% 对N进行分解N=2^m*M
[m,M]=factorize(N);
two_m=N/M;
%% 变换
if m~=0
%% 如果N可以被2整除
adr=address(m,M,two_m);
y=x(adr+1); % 蝶形运算级数l=m 时
if M~=1 % N 分解含有非2因数M时,对y中每M个数据做直接傅里叶变换for ii=1:two_m
y((ii-1)*M+1 : ii*M ) = DDFT( y((ii-1)*M+1 : ii*M ) );
end
end
%% 计算W向量
W=exp(-2*pi*i* ( 0:N/2-1 ) /N);
%% 蝶形运算
d=M;
t=N/2/d;
for l=m:-1:1
% 乘
for r=0:d-1
mul=W(r*t+1);
for ii=0:t-1
y(ii*2*d+d+1+r) = y(ii*2*d+d+1+r)*mul;
end
end
% 加
for ii=0:t-1
ind1=ii*2*d+1;
ind2=ind1+d;
for r=0:d-1
tem=y(ind1)+y(ind2);
y(ind2)=y(ind1)-y(ind2);
y(ind1)=tem;
ind1=ind1+1;
ind2=ind2+1;
end
end
d=2*d;t=t/2;
end
else
%% 如果N 不能被2整除
y=DDFT(x);
end
end
%% 内嵌函数====================================================== function y=DDFT(x)
%% 直接离散傅里叶变换
%input:
% x -- 样本数据,N维向量
%output:
% y -- N维向量
%参考文献:
% 结构动力学,克拉夫,P82
% variable define
% s -- sum,用于求和
N=length(x);
y=zeros(size(x));
for n=1:N
s=0;
for m=1:N
s=s+x(m)*exp( -i*2*pi*(m-1)*(n-1)/N );
end
y(n)=s;
end
end
function [m,M]=factorize(N)
%% 对N分解
m=0;
while true
if mod(N,2)==0
m=m+1;
N=N/2;
else
M=N;
break;
end
end
end
function adr=address(m,M,two_m)
%% 变址
% b -- 2^m * m 的矩阵,用来存储二进制数据
% ds -- 数,公差
adr=zeros(two_m,M);
b=de2bi(0:two_m-1,m);%转换为2进制注:matlab中二进制[0 1 1]=6 b=b(:,end:-1:1);% 逆序
adr(:,1)=bi2de(b);%2进制转换为10进制
if M~=1
ds=two_m;
adr=adr(:,1)*ones(1,M);
adr=adr+ds*ones(size(adr,1),1)*(0:M-1);
adr=reshape(adr',1,[]);
end
end
联系方式:matrixsuper@。