1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
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
|
