function y=extr_EKT(x,f0,Fe,m,M,K,u)
% 
%extr_EKT(x,f0,Fe,m,M,K,u) applies the EKT method to x:
% 
% Reference: J. S. Armand Eyebe Fouda,  Martin Kom and Alain Tiedeu, 
% The EKT method for extraction of periodic signals from noise, Digital 
% Signal Processing 16 (2006) 343-357.
% 
% Input:
%       x: signal to be filtered
%       f0: fundamental frequency of x
%       Fe: sampling frequency
%       m: number of desired harmonics in the output result
%       M: number of observation period for noise reduction
%       K: number of desired samples in the output signal
%       u: mean value of the reference signal (u=0 or u=1)
% 
%Output:
% 
%       y: noise reduced signal using the EKT method.
% 
% -----------------------------------------------------------------------
%                COPYRIGHT (2016): J. S. Armand EYEBE FOUDA
% -----------------------------------------------------------------------

ref=reference(Fe,f0,m,M,u);
N=length(ref);
if K>N
    K=N;
end
nx=length(x);
if nx>=N
    x1=x(1:N);
else
    x1=zeros(1,N);
    x1(1:nx)=x;
end
uu=size(x1);
y=zeros(1,K);
for k=1:K
    sr=[ref(N-k+2:N),ref(1:N-k+1)];
    if uu(1)>uu(2)
        size(sr)
        y(k)=sr*x1;
    else
        y(k)=sr*x1';
    end
end

%-----------------------------------------------------------------------

function y=reference(Fe,f0,m,M,u)
% 
%Fe est la fréquence d'échantillonnage
%f0 est la fréquence fondamentale du signal à générer,m son nombre
%d'harmoniques et M le nombre de périodes d'observation
%La fonction reference renvoie en sortie un signal de référence périodique
%utilisable par la méthode EKT.

f0=Fe/((round(M*Fe/f0))/M);
t=0:1/Fe:M/f0-1/Fe;
y=Et(m,f0,t,u)/(2*m+1);

% ----------------------------------------------------------------

function st=Et(m,f1,t,u)
b=pi*f1*t;
%u=0 si la valeur moyenne est nulle et u=1 dans le cas contraire
st=((2*m+1)*(abs(cos(b))==1))+sin((2*m+1)*b)./(sin(b)+(abs(cos(b))==1))-(u==0);
