分数阶傅立叶变换程序汇总（转载）

转载链接http://forum.vibunion.com/forum.php?mod=viewthread&action=printable&tid=41561

声振论坛
标题: 分数阶傅立叶变换程序汇总（个人收集自网上） [打印本页]
作者: 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 = [zeros(N-1,1) ; interp(f) ; zeros(N-1,1)];
% 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([x(:);y(:)])-1;
P = 2^nextpow2(N);
z = ifft( fft(x,P) .* fft(y,P));
z = z(1:N);

程序来自于http://www.cs.kuleuven.ac.be/cwi … oftware/FRFT/frft.m
作者: simon21 时间: 2007-4-26 19:51
分数阶fourior变换的源码，主要用来处理线性调频信号

%DISCRETE FRACTIONAL FOURIER TRANSFORM MATRIX GENERATOR
%by Cagatay Candan <[email]candan@ieee.org[/email]>, 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),
    [Evec,Eval]=dis_s(N,ord);
    ordp=ord;
end;

even=~rem(N,2);
F=Evec*diag(exp(-j*pi/2*a*([0:N-2 N-1+even])))*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=[dum(N); dum(1:N-1)];
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=[1 -2 1];s=0;
for k=1:ord,
    s=(-1)^(k-1)*prod(1:(k-1))^2/prod(1:2*k)*2*[0 dum(k+2:2*k+1) zeros(1,N-1-2*k) dum(1:k)]+s;
    dum=conv(dum,[1 -2 1]);
end;
S=cconvm(s)+diag(real(fft(s)));

function [Evec,Eval]=dis_s(N,ord)
%function [Evec,Eval]=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*([0:N-1])))+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);

[vc,ec]=eig(C2);[vs,es]=eig(S2);
qvc=[vc ;zeros(ceil(N/2-1),floor(N/2+1))];
SC2=P*qvc;    %Even Eigenvector of S

qvs=[zeros(floor(N/2+1),ceil(N/2-1));vs];
SS2=P*qvs;    %Odd Eigenvector of S

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

[d,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(:,[0:2:N]+1)=SC2;
    S2C2(:,[1:2:N]+1)=SS2;
    S2C2(:,N)=[];
else
    S2C2=zeros(N,N);
    S2C2(:,[0:2:N]+1)=SC2;
    S2C2(:,[1:2:N-1]+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*[0:N-1]).'.*(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 = [1 -2 1]; d_p = 1; s = 0; st = zeros(1,N);
for k = 1:p/2,
    d_p = conv(d2,d_p);
    st([N-k+1:N,1:k+1]) = 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 = [s(N:-1:2).';s(:)];
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);
[ve,ee]=eig(Ev);                 [vo,eo]=eig(Od);
%malab eig returns sorted eigenvalues
%if different routine gives unsorted eigvals, then sort first
%[d,inde] = sort(diag(ee));      [d,indo] = 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 = [1:r+1;r+2:2*r+2]; ind = ind(:);
if (even), ind([N,N+2])=[]; 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*([0:N-2 N-1+even])).'.*(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 = [1 -2 1]; d_p = 1; s = 0; st = zeros(1,N);
for k = 1:p/2,
    d_p = conv(d2,d_p);
    st([N-k+1:N,1:k+1]) = d_p; st(1) = 0;
    temp = [1:k;1:k]; temp=temp(:)'./[1:2*k];
    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 = [s(N:-1:2).';s(:)];
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);
[ve,ee]=eig(Ev);                 [vo,eo]=eig(Od);
%malab eig returns sorted eigenvalues
%if different routine gives unsorted eigvals, then sort first
%[d,inde] = sort(diag(ee));      [d,indo] = 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 = [1:r+1;r+2:2*r+2]; ind = ind(:);
if (even), ind([N,N+2])=[]; 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([fy,fys,fyss]));
    title(['a = ',num2str(a)]);
    subplot(312);plot(x,imag([fy,fys,fyss]));
    subplot(313);plot(x,abs([fy,fys,fyss]));
    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, 2001

function 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*([0:N-1])).'.*(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 = [1 -2 1]; d_p = 1; s = 0; st = zeros(1,N1);
for k = 1:p/2,
    d_p = conv(d2,d_p);
    st([N1-k+1:N1,1:k+1]) = d_p; st(1) = 0;
    temp = [1:k;1:k]; temp = temp(:)'./[1:2*k];
    s = s + (-1)^(k-1)*prod(temp)*2*st;
end;
H = toeplitz(s(:),s)+diag(real(fft(s)));
% Construct transformation matrix V
V = [zeros(N-2,1),eye(N-2),zeros(N-2,1),flipud(eye(N-2))]/sqrt(2);
V = [1,zeros(1,N1-1);V;zeros(1,N-1),1,zeros(1,N-2)];
% Compute eigenvectors
Ev = V*H*V';
[ve,ee]=eig(Ev);
% malab eig returns sorted eigenvalues
% if different routine gives unsorted eigvals, then sort first
% [d,inde] = 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, 2001

function 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*([0:N-1])).'.*(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 = [1 -2 1]; d_p = 1; s = 0; st = zeros(1,N1);
for k = 1:p/2,
    d_p = conv(d2,d_p);
    st([N1-k+1:N1,1:k+1]) = d_p; st(1) = 0;
    temp = [1:k;1:k]; temp = temp(:)'./[1:2*k];
    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 = [zeros(N,1),flipud(eye(N)),zeros(N,1),-eye(N)]/sqrt(2);
% Compute eigenvectors
Od = V*H*V';
[vo,eo]=eig(Od);
% malab eig returns sorted eigenvalues
% if different routine gives unsorted eigvals, then sort first
% [d,inde] = 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=[0:35];
y1(1:1)=1;
y1(2:36)=zeros(1,35);
subplot(323),
plot(x,y1)
axis([0,35,0,1])
title('half delta')

yc=Disfrct(y1,5/6);
subplot(324),
plot(x,real(yc),'b',x,imag(yc),'r');
axis([0,35,0-.2,0.2]);
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[signal,a_new,fact]=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;