Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

Diff of /fourier/gga.m [971186] .. [b59441] Maximize Restore

  Switch to side-by-side view

--- a/fourier/gga.m
+++ b/fourier/gga.m
@@ -1,25 +1,29 @@
-function c = gga(f,indvec,dim)
+function c = gga(f,fvec,fs,dim)
 %GGA Generalized Goertzel algorithm
 %   Usage:  c = gga(x,indvec)
 %
 %   Input parameters:
 %         x      : Input data.
-%         indvec : Indices to calculate. 
+%         fvec   : Indices to calculate. 
 %         fs     : Sampling frequency.
 %
 %   Output parameters:
 %         c      : Coefficient vector.
 %
-%   `c=gga(f,indvec)` computes the discrete-time fourier transform DTFT of
-%   *f* at 'indices' contained in `indvec`, using the generalized second-order
-%   Goertzel algorithm. Thanks to the generalization, the 'indices' can be 
-%   non-integer valued in the range 0 to *Ls-1*, where *Ls* is the length of
-%   the first non-singleton dimension of *f*. Index 0 corresponds to the
-%   DC component and integers in `indvec` result in the classical DFT
-%   coefficients. If `indvec` is empty or ommited, `indvec` is assumed to be
-%   `0:Ls-1`.
+%   `c=gga(f,fvec)` computes the discrete-time fourier transform DTFT of
+%   *f* at frequencies in `fvec` as $c(k)=F(2*pi*fvec(k))$ where
+%   $F=DTFT(f)$, $k=1,\dots\K$ and `K=length(fvec)` using the generalized
+%   second-order Goertzel algorithm. Thanks to the generalization, values
+%   in `fvec` can be arbitrary numbers in range $0-1$ and not restricted to
+%   $l/Ls$, $l=0,\dots\Ls-1$ (usual DFT samples) as the original Goertzel 
+%   algorithm is. *Ls* is the length of the first non-singleton dimension
+%   of *f*. If `indvec` is empty or ommited, `indvec` is assumed to be
+%   `(0:Ls-1)/L` and results in the same output as `fft`.
 %
-%   `c=gga(f,indvec,dim)` computes the DTFT samples along the dimension `dim`.
+%   `c=gga(f,fvec,fs)` computes the same with frequencies relative to *fs*.
+%
+%   The input *f* is processed along the first non-singleton dimension or
+%   along dimension *dim* if specified.
 %
 %   **Remark:**
 %   Besides the generalization the algorithm is also shortened by one
@@ -31,39 +35,46 @@
 %   Calculating DTFT samples of interest:::
 % 
 %     % Generate input signal
-%     k = 0:2^10-1;
-%     f = 5*sin(2*pi*k*0.05 + pi/4) + 2*sin(2*pi*k*0.1031 - pi/3);
-%
-%     % Non-integer indices of interest
-%     kgga = 102.9:0.05:109.1;
-%     % For the purposes of plot, remove the integer-valued elements
-%     kgga = setdiff(kgga,k);
-%
+%     fs = 8000;
+%     L = 2^10;
+%     k = (0:L-1).';
+%     freq = [400,510,620,680,825];
+%     phase = [pi/4,-pi/4,-pi/8,pi/4,-pi/3];
+%     amp = [5,3,4,1,2];
+%     f = arrayfun(@(a,f,p) a*sin(2*pi*k*f/fs+p),...
+%                  amp,freq,phase,'UniformOutput',0);
+%     f = sum(cell2mat(f),2);
+% 
 %     % This is equal to fft(f)
-%     ck = gga(f,k);
-%
+%     ck = gga(f);
+% 
 %     %GGA to FFT error:
 %     norm(ck-fft(f))
-%
-%     % DTFT samples just for non-integer indices
-%     ckgga = gga(f,kgga);
-%
+% 
+%     % DTFT samples at 400,510,620,680,825 Hz
+%     ckchzt = gga(f,freq,fs);
+% 
 %     % Plot modulus of coefficients
 %     figure(1);
 %     hold on;
-%     stem(k,abs(ck),'k');
-%     stem(kgga,abs(ckgga),'r:');
-%     limX = [102.9 109.1];
-%     set(gca,'XLim',limX);
-%     set(gca,'YLim',[0 1065]);
-%
+%     stem(k/L*fs,abs(ck),'k');
+%     stem(freq,abs(ckchzt),'r:');
+%     set(gca,'XLim',[freq(1)-50,freq(end)+50]);
+%     set(gca,'YLim',[0 3*L]);
+%     xlabel('f[Hz]');
+%     ylabel('|ck|');
+% 
 %     % Plot phase of coefficients
 %     figure(2);
 %     hold on;
-%     stem(k,angle(ck),'k');
-%     stem(kgga,angle(ckgga),'r:');
-%     set(gca,'XLim',limX);
+%     stem(k/L*fs,angle(ck),'k');
+%     stem(freq,angle(ckchzt),'r:');
+%     set(gca,'XLim',[freq(1)-50,freq(end)+50]);
 %     set(gca,'YLim',[-pi pi]);
+%     xlabel('f[Hz]');
+%     ylabel('angle(ck)');
+%
+%   See also: chirpzt
 %
 %   References: syra2012goertzel
        
@@ -80,23 +91,27 @@
     error('%s: X must be a nonempty vector or a matrix.',upper(mfilename))
 end
 
-if nargin<3
+if nargin<4
   dim=[];  
+end;
+
+if nargin<3 || isempty(fs)
+  fs=1;  
 end;
 
 [f,~,Ls,~,dim,permutedsize,order]=assert_sigreshape_pre(f,[],dim,'GGA');
 
-if nargin > 1 && ~isempty(indvec)
-   if ~isreal(indvec) || ~isvector(indvec)
+if nargin > 1 && ~isempty(fvec)
+   if ~isreal(fvec) || ~isvector(fvec)
       error('%s: INDVEC must be a real vector.',upper(mfilename))
    end
 else
-   indvec = 0:Ls-1;
+   fvec = (0:Ls-1)/Ls;
 end
 
-c = comp_gga(f,indvec);
+c = comp_gga(f,fvec/fs*Ls);
 
-permutedsize(1)=numel(indvec);
+permutedsize(1)=numel(fvec);
 
 c=assert_sigreshape_post(c,dim,permutedsize,order);