function y = estFreq(x,freqMuestreo)
% ESTFREQ Estimates the periodicity of a signal whose frequency is unknown.
%
% [Y] = ESTFREQ(X,FREQMUESTREO) calculates the frequency of the signal
% (Y) using the input signal (X) and the sampling rate in Hz's
% (FREQMUESTREO). The function uses autocorrelation twice.
% Copyright (C) Iván López Espejo 2010
% Version: Id: estFreq.m , v 1.0 2010
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You can obtain a copy of the GNU General Public License from
% ftp://prep.ai.mit.edu/pub/gnu/COPYING-2.0 or by writing to
% Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimation of the signal autocorrelation.
acr = zeros(1,length(x));
for k = 1:1:length(x)+1
for n = 1:1:length(x)-k
acr(k) = acr(k) + x(n)*x(k+n);
end
end
% Finally, for more precision, we study the autocorrelation of the signal
% autocorrelation.
acr2 = zeros(1,length(acr));
for k = 1:1:length(acr)+1
for n = 1:1:length(acr)-k
acr2(k) = acr2(k) + acr(n)*acr(k+n);
end
end
% Calculation of the period using the resulting function above.
OK = 0;
[max1,pos1] = max(acr2);
desp = pos1+1;
while OK ~= 1
[max2,pos2] = max(acr2(desp:length(acr2)));
pos2 = pos2 + desp;
if pos2 > 2 && pos2+1 < length(acr2)
dif1 = acr2(pos2-1) - acr2(pos2-2);
dif2 = acr2(pos2) - acr2(pos2-1);
if dif1*dif2 < 0
dif = abs(pos2-pos1);
OK = 1;
else
desp = pos2+1;
end
end
end
% Return the result.
y = freqMuestreo/dif;