快速傅里叶变换_蝶形运算_按时间抽取基2-fft算法_MATLAB代码
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
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 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@