% INFINITE TAYLOR APPROXIMATION TO DFT EIGENVECTORS
% by Ahmet Serbes <ahmet.serbes@gmail.com>, Nov. 2oo9, Istanbul
% Ahmet Serbes
%
% This code may be used for scientific and educational purposes 
% provided that credit is given to the publications below:
%        Serbes, A., & Durak-Ata, L. (2011). Efficient computation of DFT 
%        commuting matrices by a closed-form infinite order approximation 
%        to the second differentiation matrix. SIGNAL PROCESSING, 91(3), 
%        pp. 582-589.
%
%
% This Matlab function generates the discrete fractional
% Fourier transform matrix originally described in:
% Ahmet Serbes, Lutfiye Durak
%
% function M=gen_DInf(N, a) returns the NxN infinite Taylor series approxima-
% tion to the DFT commuting matrix of order a
%
% function [FrFTMtx] = gen_DInf(N,ord)
%       N   : Dimension of the FrFT matrix   
%      ord  : Fractional Fourier transform order.
%
%
% function [FrFTMtx v] = gen_DInf(N,ord)
%       N   : Dimension of the FrFT matrix   
%      ord  : Fractional Fourier transform order.
%        v  : Hermite-Gaussian-like eigenvectors.

function [FrFTMtx v] = FrFT_mtx_inf(N, ord)

% We assume that if the number of arguments at the input is 1, ord = 0.5
if(nargin == 1)
    ord = 0.5;
end

% Is N integer
if mod(N, 1) || N <= 0
    error('Dimension of the FrFT matrix has to be a positive integer');
end

% Is N even?
if(mod(N,2)==0)
    even=1;
else
    even=0;
end

% Generate the second differentiation matrix eigenvalues: df
d     =  [-2 1 zeros(1,N-3) 1];
df    =  real(fft(d));

% Infinite order differentiation matrix eigenvalues: dfinf
dfinf =  - acos((df+2)/2).^2;

% Infinite order differentiation matrix: D
u     = dft_matrix(N); % the DFT matrix...
D     = u*diag(dfinf)*u^-1;

% Infinite order commuting matrix: M
M     = D+diag(dfinf);

% Find the eigenvalues and eigenvectors of the infinite order commuting
% matrix
[v d] = eig(-M);

% Sort the eigenvalues and the corresponding eigenvectors from highest-to-
% lowest using the function sortem
[v d] = sortem(v,d);

% Now sort the HGL eigenvectors from lowest order to highest order
v=fliplr(v);

% Ignore the imaginary part. Imaginary part comes from computation errors
v = real(v);

% The eigenvalues of the FrFT
ev = diag(exp(-1i*ord*pi/2*([0:N-2 N-1+even])));

% Infinite-order approximation to the FrFT matrix.
FrFTMtx = v * ev * v';

end

function [W]=dft_matrix(N)
% W = dft_matrix(N)
% returns the normalized NxN DFT matrix 

for n=1:N
    for m=1:N
        W(n,m) = exp(j*2*(pi/N)*(n-1)*(m-1));
    end
end

% Normalization
W=W/sqrt(N);
end

function [NV,ND] = sortem(V,D)
%[V,D] = SORTEM(V,D)
%    Assumes the columns of V are vectors to be sorted along with the
%    diagonal elements of D.
%
% Matthew Dailey, 1998

% Check arguments

if nargin ~= 2
 error('Must specify vector matrix and diag value matrix')
end;

dvec = diag(D);
NV = zeros(size(V));
[dvec,index_dv] = sort(dvec);
index_dv = flipud(index_dv);
for i = 1:size(D,1)
  ND(i,i) = D(index_dv(i),index_dv(i));
  NV(:,i) = V(:,index_dv(i));
end;

end