Home > classes > @ao > computeDFT.m

computeDFT

PURPOSE ^

COMPUTEDFT Computes DFT using FFT or Goertzel

SYNOPSIS ^

function [Xx,f] = computeDFT(xin,nfft,varargin)

DESCRIPTION ^

 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 $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % COMPUTEDFT Computes DFT using FFT or Goertzel
0002 %   This function is used to calculate the DFT of a signal using the FFT
0003 %   or the Goertzel algorithm.
0004 %
0005 %   [XX,F] = COMPUTEDFT(XIN,NFFT) where NFFT is a scalar and computes the
0006 %   DFT XX using FFT. F is the frequency points at which the XX is
0007 %   computed and is of length NFFT.
0008 %
0009 %   [XX,F] = COMPUTEDFT(XIN,F) where F is a vector with atleast two
0010 %   elements computes the DFT XX using the Goertzel algorithm.
0011 %
0012 %   [XX,F] = COMPUTEDFT(...,Fs) returns the frequency vector F (in hz)
0013 %   where Fs is the sampling frequency
0014 %
0015 %   Inputs:
0016 %   XIN is the input signal
0017 %   NFFT if a scalar corresponds to the number of FFT points used to
0018 %   calculate the DFT using FFT.
0019 %   NFFT if a vector corresponds to the frequency points at which the DFT
0020 %   is calculated using goertzel.
0021 %   FS is the sampling frequency
0022 %
0023 % A direct copy of MATLAB's function for LTPDA
0024 %
0025 % M Hewitson 08-05-08
0026 %
0027 % $Id: computeDFT.m,v 1.2 2008/08/01 13:19:42 ingo Exp $
0028 %
0029 
0030 % Copyright 2006 The MathWorks, Inc.
0031 
0032 % [1] Oppenheim, A.V., and R.W. Schafer, Discrete-Time Signal Processing,
0033 % Prentice-Hall, Englewood Cliffs, NJ, 1989, pp. 713-718.
0034 % [2] Mitra, S. K., Digital Signal Processing. A Computer-Based Approach.
0035 % 2nd Ed. McGraw-Hill, N.Y., 2001.
0036 
0037 function [Xx,f] = computeDFT(xin,nfft,varargin)
0038 
0039   error(nargchk(2,3,nargin,'struct'));
0040   if nargin > 2,
0041     Fs = varargin{1};
0042   else
0043     Fs = 2*pi;
0044   end
0045 
0046   nx = size(xin,1);
0047 
0048   if length(nfft) > 1,
0049     isfreqVector = true;
0050   else
0051     isfreqVector = false;
0052   end
0053 
0054   if ~isfreqVector,
0055     [Xx,f] = computeDFTviaFFT(xin,nx,nfft,Fs);
0056   else
0057     [Xx,f] = computeDFTviaGoertzel(xin,nfft,Fs);
0058   end
0059 
0060 end
0061 
0062 %-------------------------------------------------------------------------
0063 function [Xx,f] = computeDFTviaFFT(xin,nx,nfft,Fs)
0064   % Use FFT to compute raw STFT and return the F vector.
0065 
0066   % Handle the case where NFFT is less than the segment length, i.e., "wrap"
0067   % the data as appropriate.
0068   xin_ncol = size(xin,2);
0069   xw = zeros(nfft,xin_ncol);
0070   if nx > nfft,
0071     for j = 1:xin_ncol,
0072       xw(:,j) = datawrap(xin(:,j),nfft);
0073     end
0074   else
0075     xw = xin;
0076   end
0077 
0078   Xx = fft(xw,nfft);
0079   f = psdfreqvec('npts',nfft,'Fs',Fs);
0080 end
0081 
0082 %--------------------------------------------------------------------------
0083 function [Xx,f] = computeDFTviaGoertzel(xin,freqvec,Fs)
0084   % Use Goertzel to compute raw DFT and return the F vector.
0085 
0086   f = freqvec(:);
0087   f = mod(f,Fs); % 0 <= f < = Fs
0088   nfld = floor(freqvec(:)/Fs);
0089   xm = size(xin,1); % NFFT
0090 
0091   % Indices used by the Goertzel function (see equation 11.1 pg. 755 of [2])
0092   fscaled = f/Fs*xm+1;
0093   k = round(fscaled);
0094 
0095   % shift for each frequency from default xm length grid
0096   deltak = fscaled-k;
0097 
0098   tempk = k;
0099   % If k > xm, fold over to the 1st bin
0100   k(tempk > xm) = 1;
0101   nfld = nfld + (tempk > xm); % Make nfld reflect fold in k because of round
0102 
0103   n = (0:xm-1)';
0104   Xx = zeros(size(k,1),size(xin,2));
0105   for kindex = 1:length(k)
0106     % We need to evaluate the DFT at the requested frequency instead of a
0107     % neighboring frequency that lies on the grid obtained with xm number
0108     % of points in the 0 to fs range. We do that by giving a complex phase
0109     % to xin equal to the offset from the frequency to its nearest neighbor
0110     % on the grid. This phase translates into a shift in the DFT by the
0111     % same amount. The Xx(k) then is the DFT at (k+deltak).
0112 
0113     % apply kernal to xin so as to evaluate DFT at k+deltak)
0114     kernel = exp(-j*2*pi*deltak(kindex)/xm*n);
0115     xin_phaseshifted = xin.*repmat(kernel,1,size(xin,2));
0116 
0117     Xx(kindex,:) = goertzel(xin_phaseshifted,k(kindex));
0118   end
0119 
0120   % DFT computed at exactly the frequencies it was requested for
0121   f = freqvec(:);
0122 end
0123 
0124

Generated on Mon 08-Sep-2008 13:18:47 by m2html © 2003