zak.m

linear time－frequency toolbox

function c=zak(f,a);
%ZAK  Zak transform
%   Usage:  c=zak(f,a);
%
%   ZAK(f,a) computes the Zak transform of f with parameter a.
%   The coefficients are arranged in an a x L/a matrix, where L is the
%   length of f.
%
%   If f is a matrix, then the transformation is applied to each column.
%   This is then indexed by the third dimension of the output.
%
%
%   Assume that c=ZAK(f,a), where f is a column vector of length L and
%   N=L/a. Then the following holds for m=0,...,a-1 and n=0,...,N-1
% 
%                          N-1  
%     c(m+1,n+1)=1/sqrt(N)*sum f(m-k*a+1)*exp(2*pi*i*n*k/N)
%                          k=0
%
%   SEE ALSO:  IZAK
%
%   REFERENCES:
%     H. Bölcskei and F. Hlawatsch. Discrete Zak transforms, polyphase transforms,
%     and applications. IEEE Trans. Signal Process., 45(4):851-866, april 1997.
%     
%     A. J. E. M. Janssen. Duality and biorthogonality for discrete-time
%     Weyl-Heisenberg frames. Unclassified report, Philips Electronics, 002/94.

error(nargchk(2,2,nargin));

if (prod(size(a))~=1 || ~isnumeric(a))
  error([callfun,': a must be a scalar']);
end;

if rem(a,1)~=0
  error([callfun,': a must be an integer']);
end;


if size(f,2)>1 && size(f,1)==1
  % f was a row vector.
  f=f(:);
end;

L=size(f,1);
W=size(f,2);
N=L/a;

if rem(N,1)~=0
  error('The parameter for ZAK must divide the length of the signal.');
end;

c=zeros(a,N,W);

for ii=1:W
  % Compute it, it can be done in one line!
  % We use a normalized DFT, as this gives the correct normalization
  % of the Zak transform.
  c(:,:,ii)=dft(reshape(f(:,ii),a,N),[],2);
end;