function X = gdft(x,a,b)
% GDFT Generalized discrete fourier transform
%   X = gdft(x,a,b)
%
%   The case of a=b=0 corresponds to dft
%   The case of a=1/2, b=0 corresponds to ofdft (aka odft)
%
%   x: input signal, follows fft interface of frame per column
%   a: shift in frequency (usually 0 or 1/2)
%   b: shift in time (usually 0 or 1/2)
%   X: gdft of x
%
%   Notes: This is the normalized gdft unlike the fft-ifft definition
%          in Matlab. This means that Parsevals theorem doesn't need
%          division by N

% ------- gdft.m -------------------------------------------
% Marios Athineos, marios@ee.columbia.edu
% http://www.ee.columbia.edu/~marios/
% Copyright (c) 2004-2005 by Columbia University.
% All rights reserved.
% ----------------------------------------------------------

% The default is the Odd Frequency DFT
if nargin < 2; a = 1/2; end
if nargin < 3; b = 0;   end

% Argument check
if a<0 | a>1; error('Shift in frequency must be between 0 and 1 (inclusive)'); end
if b<0 | b>1; error('Shift in time must be between 0 and 1 (inclusive)');      end

% Get size of input signal
[N,fnum] = size(x);

% Pre weights (in time)
nw = exp(-j*2*pi*a*[0:N-1]/N);
nw = diag(sparse(nw));

% Post weights (in frequency)
kw = exp(-j*2*pi*([0:N-1]+a)*b/N);
kw = diag(sparse(kw));

% So the gdft is simply (normalized)
X = full(kw*fft(full(nw*x)))/sqrt(N);