simon21 发表于 2007-4-26 19:46

分数阶傅立叶变换程序汇总(个人收集自网上)

都是从网上收集来的,由于时间比较久,处处都忘记了,如果是谁的原创请和我联系,我在帖子中标出来的

第二楼:the fast fractional Fourier transform
程序参考下文中的算法
H.M. Ozaktas, M.A. Kutay, and G. Bozdagi.
Digital computation of the fractional Fourier transform.
IEEE Trans. Sig. Proc., 44:2141--2150, 1996.

第三楼:分数阶fourior变换的源码,主要用来处理线性调频信号

第四楼:离散分数阶傅立叶变换,算法参考
C. Candan, M.A. Kutay, and H.M. Ozaktas.
The discrete Fractional Fourier Transform.
IEEE Trans. Sig. Proc., 48:1329--1337, 2000

第五楼:the discrete fractional Fourier transform,算法参考:
S.-C. Pei, M.-H. Yeh, and C.-C. Tseng.
Discrete fractional Fourier-transform based on orthogonal projections.
IEEE Trans. Sig. Proc., 47(5):1335--1348, 1999.

第六楼:a demonstration programfor the previous three routines

第七楼: the discrete fractional Cosine transform,算法参考:
S.-C. Pei, M.-H. Yeh.
Discrete fractional Cosine and Sine transforms.
IEEE Trans. Sig. Proc., 49(6):1198--1207, 2001

第八楼:the discrete fractional Sine transform,算法参考:
S.-C. Pei, M.-H. Yeh.
Discrete fractional Cosine and Sine transforms.
IEEE Trans. Sig. Proc., 49(6):1198--1207, 2001

第九楼:test for Disfrct.m,算法参考:
S.-C. Pei, M.-H. Yeh.
Discrete fractional Cosine and Sine transforms.
IEEE Trans. Sig. Proc., 49(6):1198--1207, 2001

第十楼:the rescaling preprocessing for the frft routine,详细情况参考
H.M. Ozaktas, M.A. Kutay, and G. Bozdagi.
Digital computation of the fractional Fourier transform.
IEEE Trans. Sig. Proc., 44:2141--2150, 1996.

[ 本帖最后由 simon21 于 2007-4-26 20:14 编辑 ]

simon21 发表于 2007-4-26 19:48

程序参考下文中的算法
H.M. Ozaktas, M.A. Kutay, and G. Bozdagi.
Digital computation of the fractional Fourier transform.
IEEE Trans. Sig. Proc., 44:2141--2150, 1996.

function Faf = frft(f, a)
% The fast Fractional Fourier Transform
% input: f = samples of the signal
%      a = fractional power
% output: Faf = fast Fractional Fourier transform
error(nargchk(2, 2, nargin));
f = f(:);
N = length(f);
shft = rem((0:N-1)+fix(N/2),N)+1;
sN = sqrt(N);
a = mod(a,4);
% do special cases
if (a==0), Faf = f; return; end;
if (a==2), Faf = flipud(f); return; end;
if (a==1), Faf(shft,1) = fft(f(shft))/sN; return; end
if (a==3), Faf(shft,1) = ifft(f(shft))*sN; return; end
% reduce to interval 0.5 < a < 1.5
if (a>2.0), a = a-2; f = flipud(f); end
if (a>1.5), a = a-1; f(shft,1) = fft(f(shft))/sN; end
if (a<0.5), a = a+1; f(shft,1) = ifft(f(shft))*sN; end
% the general case for 0.5 < a < 1.5
alpha = a*pi/2;
tana2 = tan(alpha/2);
sina = sin(alpha);
f = ;
% chirp premultiplication
chrp = exp(-i*pi/N*tana2/4*(-2*N+2:2*N-2)'.^2);
f = chrp.*f;
% chirp convolution
c = pi/N/sina/4;
Faf = fconv(exp(i*c*(-(4*N-4):4*N-4)'.^2),f);
Faf = Faf(4*N-3:8*N-7)*sqrt(c/pi);
% chirp post multiplication
Faf = chrp.*Faf;
% normalizing constant
Faf = exp(-i*(1-a)*pi/4)*Faf(N:2:end-N+1);

function xint=interp(x)
% sinc interpolation
N = length(x);
y = zeros(2*N-1,1);
y(1:2:2*N-1) = x;
xint = fconv(y(1:2*N-1), sinc([-(2*N-3):(2*N-3)]'/2));
xint = xint(2*N-2:end-2*N+3);

function z = fconv(x,y)
% convolution by fft
N = length()-1;
P = 2^nextpow2(N);
z = ifft( fft(x,P) .* fft(y,P));
z = z(1:N);

程序来自于http://www.cs.kuleuven.ac.be/cwis/research/nalag/research/software/FRFT/frft.m

simon21 发表于 2007-4-26 19:51

分数阶fourior变换的源码,主要用来处理线性调频信号

%DISCRETE FRACTIONAL FOURIER TRANSFORM MATRIX GENERATOR
%by Cagatay Candan <candan@ieee.org>, July 1998, Ankara
%Copyright 1998 Cagatay Candan
%This code may be used for scientific and educational purposes
%provided credit is given to the publications below:
%
%This Matlab function generates the discrete fractional
%Fourier transform matrix originally described in:
%Cagatay Candan, M. Alper Kutay, and Haldun M. Ozaktas,
%The discrete fractional Fourier Transform,
%IEEE Transactions on Signal Processing, 48:1329-1337, May 2000,
%(also in Proc ICASSP'99, pages 1713-1716, IEEE, 1999);
%and further described in:
%Haldun M. Ozaktas, Zeev Zalevsky, and M. Alper Kutay,
%The Fractional Fourier Transform with Applications in Optics and
%Signal Processing, Wiley, 2000, chapter 6.

function F=dFRT(N,a,ord)
%function F=dFRT(N,a,ord) returns the NxN discrete fractional
%Fourier transform matrix with fractional order 'a'.
%The optional argument 'ord' is the order of approximation
%of the S matrix (default 2).

%Note: This Matlab file has some subfunctions for generating S_{2k}
%      matrices, eigenvector ordering etc. These functions are not
%      visible from the Matlab workspace.

global Evec Eval ordp

if nargin==2, ord=2;end;

if (length(Evec)~=N | ordp~=ord),
    =dis_s(N,ord);
    ordp=ord;
end;

even=~rem(N,2);
F=Evec*diag(exp(-j*pi/2*a*()))*Evec';

%%%%%%

function M=cconvm(v);
%Generates circular Convm matrix

v=v(:);N=length(v);
M=zeros(N,N);dum=v;
for k=1:N,
    M(:,k)=dum;
    dum=;
end;

%%%%%%

function S=creates(N,ord)
%Creates S matrix of approximation order ord
%When ord=1, elementary S matrix is returned

ord=ord/2;
dum=;s=0;
for k=1:ord,
    s=(-1)^(k-1)*prod(1:(k-1))^2/prod(1:2*k)*2*+s;
    dum=conv(dum,);
end;
S=cconvm(s)+diag(real(fft(s)));

%%%%%%

function =dis_s(N,ord)
%function =dis_s(N)
%Returns sorted eigenvectors and eigenvalues of corresponding vectors

if nargin==1, ord=2;end;

%%Construct S Matrix
%S=diag(2*cos(2*pi/N*()))+diag(ones(1,N-1),1)+diag(ones(1,N-1),-1);
%S(1,N)=1;S(N,1)=1;
%%
S=creates(N,ord);

%%%%%%

%Construct P matrix

p=N;
r=floor(p/2);
P=zeros(p,p);

P(1,1)=1;
even=~rem(p,2);
for k=1:r-even,
    P(k+1,k+1)=1/sqrt(2);
    P(k+1,p-(k+1)+2)=1/sqrt(2);
end;

if (even), P(r+1,r+1)=1; end;

for k=r+1:p-1,
    P(k+1,k+1)=-1/sqrt(2);
    P(k+1,p-(k+1)+2)=1/sqrt(2);
end;

%%%%%%

CS=P*S*P';C2=CS(floor(1:N/2+1),floor(1:N/2+1));S2=CS(floor(N/2+2):N,floor(N/2+2):N);

=eig(C2);=eig(S2);
qvc=;
SC2=P*qvc;    %Even Eigenvector of S

qvs=;
SS2=P*qvs;    %Odd Eigenvector of S

es=diag(es);ec=diag(ec);
=sort(-ec);
SC2=SC2(:,in);
ec=ec(in);

=sort(-es);
SS2=SS2(:,in);
es=es(in);

if rem(N,2)==0,
    S2C2=zeros(N,N+1);
    SS2(:,size(SS2,2)+1)=zeros(N,1);
    S2C2(:,+1)=SC2;
    S2C2(:,+1)=SS2;
    S2C2(:,N)=[];
else
    S2C2=zeros(N,N);
    S2C2(:,+1)=SC2;
    S2C2(:,+1)=SS2;
end;

Evec=S2C2;
Eval=(-j).^[ 0:N-2 (N-1)+even];

%END

simon21 发表于 2007-4-26 19:56

离散分数阶傅立叶变换- Separate score step Fourier transforms
算法参考
C. Candan, M.A. Kutay, and H.M. Ozaktas.
The discrete Fractional Fourier Transform.
IEEE Trans. Sig. Proc., 48:1329--1337, 2000

function Fa=DFPei(f,a)

% Compute discrete fractional Fourier transform
% of order a of vector f according to Pei/Yeh/Tseng

N=length(f); f=f(:);
shft = rem((0:N-1)+fix(N/2),N)+1;

global hn_saved p_saved
if (nargin==2), p = 2; end;
p = min(max(2,p),N-1);

if (length(hn_saved) ~= N | p_saved ~= p),
    hn = make_hn(N,p);
    hn_saved = hn; p_saved = p;
else
    hn = hn_saved;
end;
Fa(shft,1)=hn*(exp(-j*pi/2*a*).'.*(hn'*f(shft)));

function hn = make_hn(N,p)
even = rem(N,2)==0;
shft = rem((0:N-1)+fix(N/2),N)+1;
% Gauss-Hermite samples
u = (-N/2:(N/2)-1)'/sqrt(N/2/pi);
ex = exp(-u.^2/2);
hn(:,1) = ex; r = norm(hn(:,1)); hn(:,1) = hn(:,1)/r;
hn(:,2)=2*u.*ex; s = norm(hn(:,2)); hn(:,2) = hn(:,2)/s;
r = s/r;
for k = 3:N+1
   hn(:,k)=2*r*u.*hn(:,k-1)-2*(k-2)*hn(:,k-2);
   s = norm(hn(:,k)); hn(:,k) = hn(:,k)/s;
   r=s/r;
end
if (even), hn(:,N)=[]; else, hn(:,N+1)=[]; end
hn(shft,:) = hn;

% eigenvectors of DFT matrix
E = make_E(N,N/2);

for k = 1:4
if even % N even
switch k
   case {1,3}
    indx = k:4:N+1;
    if (rem(N,4) ~= 0 && k==3) || (rem(N,4) == 0 && k==1),
      indx(end) = indx(end)-1;
    end
   case {2,4}
    indx = k:4:N-1;
   end
else % N odd
   indx = k:4:N;
end
hn(:,k:4:N) = orth(E(:,indx)*E(:,indx)'*hn(:,indx));
end

function E = make_E(N,p)
%Returns sorted eigenvectors and eigenvalues of corresponding vectors

%Construct matrix H, use approx order ord

d2 = ; d_p = 1; s = 0; st = zeros(1,N);
for k = 1:p/2,
    d_p = conv(d2,d_p);
    st() = d_p; st(1) = 0;
    s = s + (-1)^(k-1)*prod(1:(k-1))^2/prod(1:2*k)*2*st;      
end;
% H = circulant + diagonal
col = (0:N-1)'; row = (N:-1:1);
idx = col(:,ones(N,1)) + row(ones(N,1),:);
st = ;
H = st(idx)+diag(real(fft(s)));

%Construct transformation matrix V

r = floor(N/2);
even = ~rem(N,2);
V1 = (eye(N-1)+flipud(eye(N-1)))/sqrt(2);
V1(N-r:end,N-r:end) = -V1(N-r:end,N-r:end);
if (even), V1(r,r)=1; end
V = eye(N); V(2:N,2:N) = V1;

% Compute eigenvectors

VHV = V*H*V';
E = zeros(N);
Ev = VHV(1:r+1,1:r+1);         Od = VHV(r+2:N,r+2:N);
=eig(Ev);               =eig(Od);
%malab eig returns sorted eigenvalues
%if different routine gives unsorted eigvals, then sort first
% = sort(diag(ee));       = sort(diag(eo));
%ve = ve(:,inde');               vo = vo(:,indo');
E(1:r+1,1:r+1) = fliplr(ve);   E(r+2:N,r+2:N) = fliplr(vo);
E = V*E;
%shuffle eigenvectors
ind = ; ind = ind(:);
if (even), ind()=[]; else ind(N+1)=[]; end
E = E(:,ind');本程序来自http://www.cs.kuleuven.ac.be/cwi ... ftware/FRFT/cdpei.m

[ 本帖最后由 simon21 于 2007-4-26 19:57 编辑 ]

simon21 发表于 2007-4-26 20:00

算法参考
S.-C. Pei, M.-H. Yeh, and C.-C. Tseng.
Discrete fractional Fourier-transform based on orthogonal projections.
IEEE Trans. Sig. Proc., 47(5):1335--1348, 1999.
function y=Disfrft(f,a,p)
% Computes discrete fractional Fourier transform
% of order a of vector x
% p (optional) is order of approximation, default N/2
N = length(f); even = ~rem(N,2);
shft = rem((0:N-1)+fix(N/2),N)+1;
f = f(:);
if (nargin==2), p = N/2; end;
p = min(max(2,p),N-1);
E = dFRFT(N,p);
y(shft,1)=E*(exp(-j*pi/2*a*()).'.*(E'*f(shft)));

function E=dFRFT(N,p)
%function E=dFRFT(N,a,p) returns the NxN eigenvectors of the
%Fourier transform matrix
%The optional argument p is the order of approximation
global E_saved p_saved
if (length(E_saved) ~= N | p_saved ~= p),
    E = make_E(N,p);
    E_saved = E; p_saved = p;
else
    E = E_saved;
end;

function E = make_E(N,p)
%Returns sorted eigenvectors and eigenvalues of corresponding vectors

%Construct matrix H, use approx order ord
d2 = ; d_p = 1; s = 0; st = zeros(1,N);
for k = 1:p/2,
    d_p = conv(d2,d_p);
    st() = d_p; st(1) = 0;
    temp = ; temp=temp(:)'./;
    s = s + (-1)^(k-1)*prod(temp)*2*st;      
end;
% H = circulant + diagonal
col = (0:N-1)'; row = (N:-1:1);
idx = col(:,ones(N,1)) + row(ones(N,1),:);
st = ;
H = st(idx)+diag(real(fft(s)));
%Construct transformation matrix V
r = floor(N/2);
even = ~rem(N,2);
V1 = (eye(N-1)+flipud(eye(N-1)))/sqrt(2);
V1(N-r:end,N-r:end) = -V1(N-r:end,N-r:end);
if (even), V1(r,r)=1; end
V = eye(N); V(2:N,2:N) = V1;
% Compute eigenvectors
VHV = V*H*V';
E = zeros(N);
Ev = VHV(1:r+1,1:r+1);         Od = VHV(r+2:N,r+2:N);
=eig(Ev);               =eig(Od);
%malab eig returns sorted eigenvalues
%if different routine gives unsorted eigvals, then sort first
% = sort(diag(ee));       = sort(diag(eo));
%ve = ve(:,inde');               vo = vo(:,indo');
E(1:r+1,1:r+1) = fliplr(ve);   E(r+2:N,r+2:N) = fliplr(vo);
E = V*E;
%shuffle eigenvectors
ind = ; ind = ind(:);
if (even), ind()=[]; else ind(N+1)=[]; end
E = E(:,ind');

[ 本帖最后由 simon21 于 2007-4-26 20:08 编辑 ]

simon21 发表于 2007-4-26 20:09

a demonstration programfor the previous three routines

x=0.0:0.02:2*pi; y =cos(x);
clear p_saved hn_saved E_saved
for a=0:0.05:4
    fy=Disfrft(y,a);
    fys=cdpei(y,a);
    fyss=frft(y,a);
    % blue,green,red,cyan
    figure(1)
    subplot(311);plot(x,real());
    title(['a = ',num2str(a)]);
    subplot(312);plot(x,imag());
    subplot(313);plot(x,abs());
    pause(0.7);
end

simon21 发表于 2007-4-26 20:10

the discrete fractional Cosine

算法参考:
S.-C. Pei, M.-H. Yeh.
Discrete fractional Cosine and Sine transforms.
IEEE Trans. Sig. Proc., 49(6):1198--1207, 2001function y = Disfrct(f,a,p)
% Computes discrete fractional cosine transform
% of order a of vector f
% p (optional) is order of approximation, default N/2
% S-C Pei, M-H Yeh, IEEE Tr SP 49(6)2001, pp.1198-1207
N = length(f);
shft = rem((0:N-1)+fix(N/2),N)+1;
f = f(:);
if (nargin==2), p = N/2; end;
p = min(max(2,p),N-1);
E = dFRCT(N,p);
y=E*(exp(-j*pi*a*()).'.*(E'*f));

function E=dFRCT(N,p)
%function E=dFRCT(N,p) returns the NxN eigenvectors of the
%Fourier Cosine transform matrix

global EC_saved pC_saved

if (length(EC_saved) ~= N | pC_saved ~= p),
    E = make_EC(N,p);
    EC_saved = E; pC_saved = p;
else
    E = EC_saved;
end;

function E = make_EC(N,p)
% Returns sorted eigenvectors and eigenvalues of corresponding vectors
% Construct matrix H, use approx order p
N1=2*N-2;
d2 = ; d_p = 1; s = 0; st = zeros(1,N1);
for k = 1:p/2,
    d_p = conv(d2,d_p);
    st() = d_p; st(1) = 0;
    temp = ; temp = temp(:)'./;
    s = s + (-1)^(k-1)*prod(temp)*2*st;
end;
H = toeplitz(s(:),s)+diag(real(fft(s)));
% Construct transformation matrix V
V = /sqrt(2);
V = ;
% Compute eigenvectors
Ev = V*H*V';
=eig(Ev);
% malab eig returns sorted eigenvalues
% if different routine gives unsorted eigvals, then sort first
% = sort(diag(ee));
% ve = ve(:,inde');
E = fliplr(ve);
E(end,:) = E(end,:)/sqrt(2);

[ 本帖最后由 simon21 于 2007-4-26 20:11 编辑 ]

simon21 发表于 2007-4-26 20:11

the discrete fractional Sine transform
算法参考:
S.-C. Pei, M.-H. Yeh.
Discrete fractional Cosine and Sine transforms.
IEEE Trans. Sig. Proc., 49(6):1198--1207, 2001function y = Disfrst(f,a,p)
% Computes discrete fractional sine transform
% of order a of vector f
% p (optional) is order of approximation, default N/2
% S-C Pei, M-H Yeh, IEEE Tr SP 49(6)2001, pp.1198-1207
N = length(f);
shft = rem((0:N-1)+fix(N/2),N)+1;
f = f(:);
if (nargin==2), p = N/2; end;
p = min(max(2,p),N-1);
E = dFRST(N,p);
y=E*(exp(-j*pi*a*()).'.*(E'*f));

function E=dFRST(N,p)
%function E=dFRST(N,p) returns the NxN eigenvectors of the
%Fourier Sine transform matrix

global ES_saved pS_saved

if (length(ES_saved) ~= N | pS_saved ~= p),
    E = make_ES(N,p);
    ES_saved = E; pS_saved = p;
else
    E = ES_saved;
end;

function E = make_ES(N,p)
% Returns sorted eigenvectors and eigenvalues of corresponding vectors
% Construct matrix H, use approx order p
N1=2*N+2;
d2 = ; d_p = 1; s = 0; st = zeros(1,N1);
for k = 1:p/2,
    d_p = conv(d2,d_p);
    st() = d_p; st(1) = 0;
    temp = ; temp = temp(:)'./;
    s = s + (-1)^(k-1)*prod(temp)*2*st;
end;
H = toeplitz(s(:),s)+diag(real(fft(s)));
% Construct transformation matrix V
r = N;
V = /sqrt(2);
% Compute eigenvectors
Od = V*H*V';
=eig(Od);
% malab eig returns sorted eigenvalues
% if different routine gives unsorted eigvals, then sort first
% = sort(diag(eo));
% vo = vo(:,inde');
E = flipud(fliplr(vo));

simon21 发表于 2007-4-26 20:13

test for Disfrct.m
算法参考:
S.-C. Pei, M.-H. Yeh.
Discrete fractional Cosine and Sine transforms.
IEEE Trans. Sig. Proc., 49(6):1198--1207, 2001

% Disfrft

function testfrct(x)
x=[-35:35];
y = zeros(size(x));
y(fix(length(x)/2)+1)=1;
subplot(321),
plot(x,y);
title('symmetric delta')

yf=Disfrft(y,5/6);
subplot(322),
plot(x,real(yf),'b',x,imag(yf),'r');
axis([-32,32,-0.2,0.2]);
title('FrFT for a = 5/6 of symmetric delta')

% Disfrct

x=;
y1(1:1)=1;
y1(2:36)=zeros(1,35);
subplot(323),
plot(x,y1)
axis()
title('half delta')

yc=Disfrct(y1,5/6);
subplot(324),
plot(x,real(yc),'b',x,imag(yc),'r');
axis();
title('FrCT for a = 5/6 of half delta')

% Disfrcst

x=[-35:0];
y1(1:36)=zeros(1,36);
y1(1:1)=-1;
subplot(325),
plot(x,y1)
axis([-35,0,-1,0])
title('half delta')

ys=Disfrst(y1,5/6);
subplot(326),
plot(x,real(ys),'b',x,imag(ys),'r');
axis([-35,0,-.2,0.2]);
title('FrST for a = 5/6 of half delta')

simon21 发表于 2007-4-26 20:14

the rescaling preprocessing for the frft routine

算法参考:
H.M. Ozaktas, M.A. Kutay, and G. Bozdagi.
Digital computation of the fractional Fourier transform.
IEEE Trans. Sig. Proc., 44:2141--2150, 1996.
function=rescale(signal,a,delta_x1);
% This routine calculates the parameters to transform the input for the
% fractional Fourier transform when the support of the input
% input is not N^2 with N=sqrt(length(signal))
% Parameters: signal   = the signal to be transformed
%             a      = the order of the transform
%             delta_x1 = the length of the support of the input signal.
% To compute the frft with these data use
% Output:   signal = possibly flipped signal to avoid infinite factor
%             a_new = compute the frft of order a_new
%             fact = and elementwise multiply the result with fact
N = length(signal);
delta_x2 = sqrt(N);
k = delta_x2/delta_x1;

a_new =mod(a,4);

if k == 1 || a == 0
a_new = a;
fact = 1;
return
elseif a == 2
signal = flipud(signal)
else
phi = a*pi/2;
psi = acot(cot(phi)/k^2);
c = csc(phi)/(k*csc(psi));
x = linspace(-delta_x1/2,delta_x1/2,N);
u = x;
nu = c*u;
a_new = 2*psi/pi;
A_phi = exp(-i*(pi*sign(sin(phi))/4-phi/2))/sqrt(abs(sin(phi)));
A_psi = exp(-i*(pi*sign(sin(psi))/4-psi/2))/sqrt(abs(sin(psi)));
fact = A_phi/(k*A_psi)*exp(i*pi*(cot(phi)/c^2-cot(psi))*(nu.^2))';
end;

wangyuangan 发表于 2008-4-4 17:50

谢谢

:handshake
:lol

xueer534 发表于 2008-4-19 12:37

回复 10楼 的帖子

我看不太懂这个程序,感觉对输入输出参数的理解也很模糊,请清楚的同行解释一下此程序吧。
在调用frft.m程序前,需要事先调用此函数吗?

dailiangren 发表于 2008-4-19 16:59

。。。。。。。。。。好帖,留名。

wyhlucky 发表于 2008-4-21 17:27

楼主有没有C的程序啊?急急急啊。

cxblhll 发表于 2008-4-29 15:14

晕死,看不懂,看来对信号这方面俺还是张白纸啊
页: [1] 2 3 4
查看完整版本: 分数阶傅立叶变换程序汇总(个人收集自网上)