--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/computeDFT.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,124 @@
+% COMPUTEDFT Computes DFT using FFT or Goertzel
+%   This function is used to calculate the DFT of a signal using the FFT
+%   or the Goertzel algorithm.
+%
+%   [XX,F] = COMPUTEDFT(XIN,NFFT) where NFFT is a scalar and computes the
+%   DFT XX using FFT. F is the frequency points at which the XX is
+%   computed and is of length NFFT.
+%
+%   [XX,F] = COMPUTEDFT(XIN,F) where F is a vector with atleast two
+%   elements computes the DFT XX using the Goertzel algorithm.
+%
+%   [XX,F] = COMPUTEDFT(...,Fs) returns the frequency vector F (in hz)
+%   where Fs is the sampling frequency
+%
+%   Inputs:
+%   XIN is the input signal
+%   NFFT if a scalar corresponds to the number of FFT points used to
+%   calculate the DFT using FFT.
+%   NFFT if a vector corresponds to the frequency points at which the DFT
+%   is calculated using goertzel.
+%   FS is the sampling frequency
+%
+% A direct copy of MATLAB's function for LTPDA
+%
+% M Hewitson 08-05-08
+%
+% $Id: computeDFT.m,v 1.2 2008/08/01 13:19:42 ingo Exp $
+%
+
+% Copyright 2006 The MathWorks, Inc.
+
+% [1] Oppenheim, A.V., and R.W. Schafer, Discrete-Time Signal Processing,
+% Prentice-Hall, Englewood Cliffs, NJ, 1989, pp. 713-718.
+% [2] Mitra, S. K., Digital Signal Processing. A Computer-Based Approach.
+% 2nd Ed. McGraw-Hill, N.Y., 2001.
+
+function [Xx,f] = computeDFT(xin,nfft,varargin)
+
+  error(nargchk(2,3,nargin,'struct'));
+  if nargin > 2,
+    Fs = varargin{1};
+  else
+    Fs = 2*pi;
+  end
+
+  nx = size(xin,1);
+
+  if length(nfft) > 1,
+    isfreqVector = true;
+  else
+    isfreqVector = false;
+  end
+
+  if ~isfreqVector,
+    [Xx,f] = computeDFTviaFFT(xin,nx,nfft,Fs);
+  else
+    [Xx,f] = computeDFTviaGoertzel(xin,nfft,Fs);
+  end
+
+end
+
+%-------------------------------------------------------------------------
+function [Xx,f] = computeDFTviaFFT(xin,nx,nfft,Fs)
+  % Use FFT to compute raw STFT and return the F vector.
+
+  % Handle the case where NFFT is less than the segment length, i.e., "wrap"
+  % the data as appropriate.
+  xin_ncol = size(xin,2);
+  xw = zeros(nfft,xin_ncol);
+  if nx > nfft,
+    for j = 1:xin_ncol,
+      xw(:,j) = datawrap(xin(:,j),nfft);
+    end
+  else
+    xw = xin;
+  end
+
+  Xx = fft(xw,nfft);
+  f = psdfreqvec('npts',nfft,'Fs',Fs);
+end
+
+%--------------------------------------------------------------------------
+function [Xx,f] = computeDFTviaGoertzel(xin,freqvec,Fs)
+  % Use Goertzel to compute raw DFT and return the F vector.
+
+  f = freqvec(:);
+  f = mod(f,Fs); % 0 <= f < = Fs
+  nfld = floor(freqvec(:)/Fs);
+  xm = size(xin,1); % NFFT
+
+  % Indices used by the Goertzel function (see equation 11.1 pg. 755 of [2])
+  fscaled = f/Fs*xm+1;
+  k = round(fscaled);
+
+  % shift for each frequency from default xm length grid
+  deltak = fscaled-k;
+
+  tempk = k;
+  % If k > xm, fold over to the 1st bin
+  k(tempk > xm) = 1;
+  nfld = nfld + (tempk > xm); % Make nfld reflect fold in k because of round
+
+  n = (0:xm-1)';
+  Xx = zeros(size(k,1),size(xin,2));
+  for kindex = 1:length(k)
+    % We need to evaluate the DFT at the requested frequency instead of a
+    % neighboring frequency that lies on the grid obtained with xm number
+    % of points in the 0 to fs range. We do that by giving a complex phase
+    % to xin equal to the offset from the frequency to its nearest neighbor
+    % on the grid. This phase translates into a shift in the DFT by the
+    % same amount. The Xx(k) then is the DFT at (k+deltak).
+
+    % apply kernal to xin so as to evaluate DFT at k+deltak)
+    kernel = exp(-j*2*pi*deltak(kindex)/xm*n);
+    xin_phaseshifted = xin.*repmat(kernel,1,size(xin,2));
+
+    Xx(kindex,:) = goertzel(xin_phaseshifted,k(kindex));
+  end
+
+  % DFT computed at exactly the frequencies it was requested for
+  f = freqvec(:);
+end
+
+