function      Y = goertzelmex(X,I)
%
%	Goertzel algorithm
%			1 - exp(-j*2*pi*k/N)*z^-1
%	Hk(z)=----------------------------------
%			1 - 2*cos(2*pi*k/N)z^-1 + z^-2
%
%	Peter Karlsen 23/10-2005
%

N=length(X);
K=length(I);

for k=1:K
  omega = (2.0 * pi * I(k)) / N;
  sine = sin(omega);
  cosine = cos(omega);
  coeff = 2.0 * cosine

  Q2 = 0;
  Q1 = 0;
  for n=1:N
          Q0 = X(n) + coeff * Q1 - Q2;
          Q2 = Q1;
          Q1 = Q0;
  end
  Y(k)=  Q1 - Q2 * exp(-j*omega);
end
