function [y,f] = dtftfft(x,varargin)
%dtftfft  Compute the discrete-time Fourier transform using FFT.
%
%   dtftfft(x) computes the dtft of
%   x = [x_0,...,x_{N-1}] evaluated for
%   f = (k-floor(N/2))/(2N) for k=0,...,N-1.
%   In other words, f is multiples of (+/-)1/N
%   that lie in [-1/2,1/2].
%
%   dtftfft(x,M1) computes the dtft of
%   x = [x_{M1},...,x_{M2}], where
%   M2 is implicitly defined as
%   length(x)+M1-1.
%
%   dtftfft(x,M1,N) computes the dtft of
%   x = [x_{M1},...,x_{M2},0,...,0], where
%   x is zero padded (or truncated) to length N.
%   If N=0, x is zero padded to a length
%   that is a power of 2.
%   If N<0, x is zero padded to a length that is
%   a power of 2 and also greater than or equal to
%   both the length of x and -N.

N = length(x);
if nargin==3
   N = varargin{2};
   if N==0 % If N=0, zero pad x to a power of 2.
      N = 2^ceil(log2(length(x)));
   elseif N<0 % If N<0, zero pad x to a power of 2 >= -N
      N = 2^ceil(log2(max(length(x),-N))); % and >= length(x)
   end
end
y = fft(x,N);
k = [0:N-1];
if nargin>=2
   M1 = varargin{1};
   if M1~=0
      y = exp(-j*2*pi*M1/N*k).*y;
   end
end
y = fftshift(y);
% For N=2m even, change 0,...,m,...,2m-1
% to -m,...,0,...,m-1.  For N=2m+1 odd, change
% 0,...,m-1,m,m+1,...,2m to -m,...,-1,0,1,...,m.
% Then divide by N to get frequencies in [-1/2,1/2].
f = (k-floor(N/2))/N;