📄 partial fast fourier transform fpft.mht
字号:
From: <由 Windows Internet Explorer 7 保存>
Subject:
Date: Sat, 26 Apr 2008 08:21:14 +0800
MIME-Version: 1.0
Content-Type: text/html;
charset="gb2312"
Content-Transfer-Encoding: quoted-printable
Content-Location: http://www.mathworks.com/matlabcentral/files/9200/fpft.m
X-MimeOLE: Produced By Microsoft MimeOLE V6.00.2900.3198
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META http-equiv=3DContent-Type content=3D"text/html; charset=3Dgb2312">
<META content=3D"MSHTML 6.00.6000.16640" name=3DGENERATOR></HEAD>
<BODY><PRE>function fpft=3Dfastpartialfouriertransform(x,n,m,shift);
% FPFT fast partial fourier transform
%
% FPFT(X) is equal to FFT(X)
% FPFT(X,N) is equal to FFT(X,N), X is padded with zeroes to length N
% FPFT(X,N,M) is equal to the first M points (rows) of the FFT(X,N)
% FPFT(X,N,M,SHIFT) is equal to the points (rows) SHIFT+1:SHIFT+M of the =
FFT (X,N)
% N, M and size(X,1) must be a power of 2
%
% For matrices FPFT (like FFT) operates along the 1st dimension of X.
% FPFT was based on the description in
% The Fractional Fourier Transform and Applications
% Daivd H. Bailey and Paul N. Swarztrauber
% October 19, 1995
% SIAM Review, vol. 33 no 3 (Sept. 1991), pg 389-404
% Download: http://crd.lbl.gov/~dhbailey/dhbpapers/
% input argument checking
if (nargin<1)
error('Not enough input arguments.')
end
if (ndims(x)>2)
error('X must be vector or 2-D matrix.');
end;
if (nargin<2) % just return the fft
fpft=3Dfft(x);
return;
end;
if (nargin<3 || isempty(m) || m=3D=3Dn) % just return the fft
fpft=3Dfft(x,n);
return;
end;
% vec=3D0 array, vec=3D1 column vector, vec=3D2 row vector
[lx wx]=3Dsize(x);
if wx~=3D1
if lx~=3D1
vec=3D0; % 2 dim array
else
x=3Dx.';
lx=3Dwx;
wx=3D1;
vec=3D2; % row vector
end
else % column vector
vec=3D1;=20
end;
if isempty(n)
n=3Dlx;
else
if (n<lx), % truncate x
x=3Dx(1:n,:);
lx=3Dn;
end;
end;
if (m>n)
error('M must be <=3D N.');
end;
if rem(log2(n),1) || rem(log2(m),1) || rem(log2(lx),1)
error('N, size(X,1) and M must be powers of 2.');
end
if nargin<4 || isempty(shift)
shift=3D0;
else
shift=3Dshift-round(shift/n)*n; % shift must be -n/2 <=3D shift =
<=3D n/2
end;
if m>=3Dn/2 || lx>=3Dn/4 % in this case just do a (more =
efficient) normal FFT
fpft=3Dfft(x,n);
=20
% this applies the circular shift
=20
range=3D[min(n+1,n+shift+1):min(n,n+shift+m) =
max(1,shift+1):min(shift+m,n) max(1,shift+m-n):max(0,shift+m-n) ];
=20
switch vec,
case 0
fpft=3Dfpft(range,:);
case 1
fpft=3Dfpft(range);
case 2;
fpft=3Dfpft(range).';
end;
=20
return;
end;
% now the real magic begins
% from here use the fast fractional fourier transform
persistent expf lookupn countn
maxn=3D10; % maximum number of cachable n's, increase this number if =
needed
if (isempty(countn)) % initialize persistent variables
lookupn=3Drepmat([NaN 0],maxn,1); % cache hit counter
expf=3Dzeros(maxn,n+1);
countn=3D0;
end;
indn=3Dfind(lookupn(:,1)=3D=3Dn);
if (isempty(indn))
if (countn>=3Dmaxn), % no more cache-able
[mn,indn]=3Dmin(lookupn(:,2)); % delete least used
else
countn=3Dcountn+1;
indn=3Dcountn;
end; =20
lookupn(indn,:)=3D[n 0];
expf(indn,1:n+1)=3Dexp(-i*pi/n*(0:n).^2);
else
lookupn(indn,2)=3Dlookupn(indn,2)+1; % increase usage count
end;
indwx=3Dones(1,wx);
indnm=3Dindn(indwx);
ffty=3Dfft(x.*expf(indnm,1:lx).',2*lx);
if (m>lx)
fpft=3Dzeros(wx,m);
=20
for s=3Dshift:lx:shift+m-1
indz=3Dabs([s:lx-1+s -lx+s:-1+s])+1;
ftz=3Dfft(1./expf(indn,indz).');
cv=3Difft(ftz(:,indwx).*ffty).';
=
fpft(:,s-shift+1:s-shift+lx)=3Dcv(:,1:lx).*expf(indnm,indz(1:lx));
end;
else
indz=3Dabs([shift:lx-1+shift -lx+shift:-1+shift])+1;
ftz=3Dfft(1./expf(indn,indz).');
cv=3Difft(ftz(:,indwx).*ffty).';
fpft=3Dcv(:,1:m).*expf(indnm,indz(1:m));
end;
if (vec~=3D2)
fpft=3Dfpft.';
end;
</PRE></BODY></HTML>
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -