function X_k = gfta(x, k)
% GFTA (Goertzel Fourier transform algorithm)
% Вычисление компонент спектра по алгоритму Герцеля
%
% X_k = GAFT(x, k)
%
% x   - входные данные               % N x 1
% k   - номера вычисляемых компонент % M x 1
%
% X_k - выходные данные              % M x 1
%
% Число умножений: M(N-1)
% Число сложений: 2M(N-1)
% При N > 6 эффективнее прямого вычисления
%
% 3 ноября 2004 г.
%
% Просеков О. В. (sc2@pisem.net) 
%
% Санкт-Петербургский городской семинар «Всплески (wavelets) и их приложения».
% Секция «Дискретный гармонический анализ»: http://www.math.spbu.ru/user/dmp/dha/

N = length(x);
theta = 2*pi/N * k;
X_k = Goertzel_algorithm(x, theta);

function X = Goertzel_algorithm(x, theta)
N = length(x);
M = length(theta);
c = 2*cos(theta);
g = zeros(M, 3);
g(:,1) = x(N);
for j = N-1:-1:2
    g(:,[3 2]) = g(:,[2 1]);
    g(:,1) = x(j) + c .* g(:,2) - g(:,3);
end
X = x(1) - g(:,2) + (cos(theta) - 1i * sin(theta)) .* g(:,1);