line source
+ − % 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
+ −
+ −