1 files changed, 79 insertions, 0 deletions
diff --git a/macros/goertzel.sci b/macros/goertzel.sci
new file mode 100644
index 0000000..bc0061c
--- /dev/null
+++ b/macros/goertzel.sci
@@ -0,0 +1,79 @@
+function Y = goertzel(X,INDVEC,DIM)
+//Computes DFT using the second order Goertzel Algorithm
+//Calling Sequence
+//Y = goertzel(X,INDVEC,DIM)
+//Parameters
+//X
+//A vector matrix or n-dimensional array
+//INDVEC
+//The indices at which the DFT is to be computed
+//DIM
+//The dimension along which the algorithm is to be implemented
+//Description
+//goertzel(X,INDVEC)
+//Computes the DFT of X at indices INDVEC using the second order algorithm along
+//the first non-singleton dimension. Elements of INDVEC must be positive integers
+//less than the length of the first non-singleton dimension. If INDVEC is empty
+//the DFT is computed at all indices along the first non-singleton dimension
+//goertzel(X,INDVEC,DIM)
+//Implements the algorithm along dimension DIM
+//In general goertzel is slower than fft when computing the DFT for all indices
+//along a particular dimension. However it is computationally more efficient when
+//the DFT at only a subset of indices is desired
+//Example
+//x=rand(1,5)
+//x  =
+// 
+//    0.6283918    0.8497452    0.6857310    0.8782165    0.0683740  
+//y=goertzel(x,2);
+//y  =
+// 
+//  - 0.3531539 - 0.6299881i  
+//Author
+//Ankur Mallick
+//References
+//Goertzel, G. (January 1958), "An Algorithm for the Evaluation of Finite Trigonometric Series", American Mathematical Monthly 65 (1): 34–35, doi:10.2307/2310304
+    funcprot(0);
+    if(argn(2)<3|isempty(DIM))
+        DIM=find(size(X)>1,1); //First non-singleton dimension
+    end
+    if(DIM>ndims(X))
+        error('Invalid Dimensions');
+    end
+    perm=[DIM, 1:DIM-1, DIM+1:ndims(X)];
+    X=permute(X,perm); //Makes DIM the leading dimension
+    S=size(X);
+    if(argn(2)<2|isempty(INDVEC))
+        INDVEC=1:S(1);
+    end
+    if(max(INDVEC)>S(1)|min(INDVEC)<1)
+        error('Index out of bounds');
+    elseif(or(INDVEC~=round(INDVEC)))
+        error('Indices must be integers');
+    end
+    X1=matrix(X,1,prod(S));
+    T=[type(X1), type(INDVEC), type(DIM)]; //all inputs should be of type double
+    if(~and(T==1))
+        error('Invalid data type');
+    end
+    //Implementing Goertzel algorithm
+    len=[length(INDVEC), S(2:length(S))];
+    Y=matrix(zeros(1,prod(len)),len);
+    v=ones(length(INDVEC),1);
+    w=(2*%pi*(INDVEC-1)/S(1))';
+    re=2*cos(w);
+    im=sin(w);
+    for i=1:prod(S(2:length(S)))
+        x=X(:,i)
+        sp1=0;
+        sp2=0;
+        for j=1:S(1)
+            s=x(j)*v+re.*sp1-sp2;
+            sp2=sp1;
+            sp1=s;
+        end
+        Y(:,i)=((re/2)+im*%i).*sp1-sp2;
+    end
+    iperm(perm)=1:length(perm);
+    Y=permute(Y,iperm); //Converting Y to the original shape of X
+endfunction