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 $
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