% For WABF 2006/2007, R. Vargic, 2006
% Comparison of different methods for computration of center frequency
% of various basic wavelets 
% Remarks: I propose to use precision=10
%          To see some nice graphics, please enable the plot on line 41
%          To use different algorithm, just edit the fnc declaration
%          To support more wavelets edit the  "if" condition on line 14

function MEFrq     = strfrq(wName,precision)
%function maxFrq    = strfrq(wName,precision)
%function matlabFrq = strfrq(wName,precision)

%============================================================
%reference from Matlab
matlabFrq=centfrq(wName,precision)

%============================================================
%time samples extension for periodification
extension=64;
if (strcmp(wName,'gaus1')==1)
    [tpsi,os]=wavefun(wName,precision);
else
    [tphi,tpsi,os]=wavefun(wName,precision);
end
%plot(os,tphi);
tsamples=length(os);
timeMin=os(1);
timeMax=os(tsamples);
tvz=(timeMax-timeMin)/tsamples;
fvz=1/tvz;
%we extend the samples number - we extend the frequency presision
psamples=tsamples*extension;
perFnc=zeros(1,psamples);
perFnc(1:tsamples)=tpsi;
%plot(perFnc);
%sfnc=fft(tfnc)
sfnc=fft(perFnc);
unitFreq=fvz/psamples;
frqscale=0:unitFreq:(fvz-1);
dispFactor=ceil(psamples/500);
%plot(frqscale(1:dispFactor),abs(sfnc(1:dispFactor)));

%============================================================
%search for maximum energy
[maxVal, maxPos]=max(abs(sfnc));
%interprete the result
maxFrq=unitFreq*maxPos


%============================================================
%search for the center of time-frequency window
energy=abs(sfnc).*abs(sfnc);
myProduct=0;
for idx = 1:psamples/2    
    myProduct=myProduct+idx*energy(idx);
end
%Mean Energy position
MEXpos=myProduct/sum(energy(2:psamples/2));
%mean energy asociated frequency
MEFrq=unitFreq*MEXpos