Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/computeDFT.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % COMPUTEDFT Computes DFT using FFT or Goertzel
2 % This function is used to calculate the DFT of a signal using the FFT
3 % or the Goertzel algorithm.
4 %
5 % [XX,F] = COMPUTEDFT(XIN,NFFT) where NFFT is a scalar and computes the
6 % DFT XX using FFT. F is the frequency points at which the XX is
7 % computed and is of length NFFT.
8 %
9 % [XX,F] = COMPUTEDFT(XIN,F) where F is a vector with atleast two
10 % elements computes the DFT XX using the Goertzel algorithm.
11 %
12 % [XX,F] = COMPUTEDFT(...,Fs) returns the frequency vector F (in hz)
13 % where Fs is the sampling frequency
14 %
15 % Inputs:
16 % XIN is the input signal
17 % NFFT if a scalar corresponds to the number of FFT points used to
18 % calculate the DFT using FFT.
19 % NFFT if a vector corresponds to the frequency points at which the DFT
20 % is calculated using goertzel.
21 % FS is the sampling frequency
22 %
23 % A direct copy of MATLAB's function for LTPDA
24 %
25 % M Hewitson 08-05-08
26 %
27 % $Id: computeDFT.m,v 1.2 2008/08/01 13:19:42 ingo Exp $
28 %
29
30 % Copyright 2006 The MathWorks, Inc.
31
32 % [1] Oppenheim, A.V., and R.W. Schafer, Discrete-Time Signal Processing,
33 % Prentice-Hall, Englewood Cliffs, NJ, 1989, pp. 713-718.
34 % [2] Mitra, S. K., Digital Signal Processing. A Computer-Based Approach.
35 % 2nd Ed. McGraw-Hill, N.Y., 2001.
36
37 function [Xx,f] = computeDFT(xin,nfft,varargin)
38
39 error(nargchk(2,3,nargin,'struct'));
40 if nargin > 2,
41 Fs = varargin{1};
42 else
43 Fs = 2*pi;
44 end
45
46 nx = size(xin,1);
47
48 if length(nfft) > 1,
49 isfreqVector = true;
50 else
51 isfreqVector = false;
52 end
53
54 if ~isfreqVector,
55 [Xx,f] = computeDFTviaFFT(xin,nx,nfft,Fs);
56 else
57 [Xx,f] = computeDFTviaGoertzel(xin,nfft,Fs);
58 end
59
60 end
61
62 %-------------------------------------------------------------------------
63 function [Xx,f] = computeDFTviaFFT(xin,nx,nfft,Fs)
64 % Use FFT to compute raw STFT and return the F vector.
65
66 % Handle the case where NFFT is less than the segment length, i.e., "wrap"
67 % the data as appropriate.
68 xin_ncol = size(xin,2);
69 xw = zeros(nfft,xin_ncol);
70 if nx > nfft,
71 for j = 1:xin_ncol,
72 xw(:,j) = datawrap(xin(:,j),nfft);
73 end
74 else
75 xw = xin;
76 end
77
78 Xx = fft(xw,nfft);
79 f = psdfreqvec('npts',nfft,'Fs',Fs);
80 end
81
82 %--------------------------------------------------------------------------
83 function [Xx,f] = computeDFTviaGoertzel(xin,freqvec,Fs)
84 % Use Goertzel to compute raw DFT and return the F vector.
85
86 f = freqvec(:);
87 f = mod(f,Fs); % 0 <= f < = Fs
88 nfld = floor(freqvec(:)/Fs);
89 xm = size(xin,1); % NFFT
90
91 % Indices used by the Goertzel function (see equation 11.1 pg. 755 of [2])
92 fscaled = f/Fs*xm+1;
93 k = round(fscaled);
94
95 % shift for each frequency from default xm length grid
96 deltak = fscaled-k;
97
98 tempk = k;
99 % If k > xm, fold over to the 1st bin
100 k(tempk > xm) = 1;
101 nfld = nfld + (tempk > xm); % Make nfld reflect fold in k because of round
102
103 n = (0:xm-1)';
104 Xx = zeros(size(k,1),size(xin,2));
105 for kindex = 1:length(k)
106 % We need to evaluate the DFT at the requested frequency instead of a
107 % neighboring frequency that lies on the grid obtained with xm number
108 % of points in the 0 to fs range. We do that by giving a complex phase
109 % to xin equal to the offset from the frequency to its nearest neighbor
110 % on the grid. This phase translates into a shift in the DFT by the
111 % same amount. The Xx(k) then is the DFT at (k+deltak).
112
113 % apply kernal to xin so as to evaluate DFT at k+deltak)
114 kernel = exp(-j*2*pi*deltak(kindex)/xm*n);
115 xin_phaseshifted = xin.*repmat(kernel,1,size(xin,2));
116
117 Xx(kindex,:) = goertzel(xin_phaseshifted,k(kindex));
118 end
119
120 % DFT computed at exactly the frequencies it was requested for
121 f = freqvec(:);
122 end
123
124