Instaneous_Frequency_Test.m

function [f,t,f_GZC]=Instaneous_Frequency_Test(x,fs)
% [f,t]=Instaneous_Frequency_Test(x,fs)[f,t]=Instaneous_Frequency_Test(x,fs)
% x is the signal that will be analysised.
% fs is the sampling frequency
% f is the instaneous frquency
% t is the time range of each instaneous frequency
% f_GZC is the mean frequency

%the length of x must be larger than 8
if length(x)<8
    error('the length of x must be larger than 7')
end
%clear results
f=[]; 
t=[];
%calculate sampling periode
Ts=1/fs;
%calculate local extrema and zero
[indmin indmax indzer]=extr(x);
%refresh local extrema and zero
[ newindmin newindmax newindzer ] = refresh_index( x,indmin,indmax,indzer );
indzer=newindzer;
indmin=newindmin;
indmax=newindmax;
%get the index of all indmin, indmax and indzer
%[index,I] = rerange_index( indmin,indmax,indzer );
index=sort([indmin indmax indzer]);
if length(index)<5
    error('Too few extrema and zero')
end
%prepare Begin point index and finish point index
L_index=length(index);
%begin calculate instaneous frquency
Begin_point=4;
Finish_point=L_index-4;
for Number_index=1:L_index-1
    if Number_index==1
        t=[t;index(Number_index) index(Number_index+1)];
        t_q=(index(Number_index+1)-index(Number_index))*Ts;
        t_h_2=(index(Number_index+2)-index(Number_index))*Ts;
        t_f_4=(index(Number_index+4)-index(Number_index))*Ts;
        f=[f 1/3*(1/t_q/4+1/t_h_2/2+1/t_f_4)];
    elseif Number_index==2
        t=[t;index(Number_index) index(Number_index+1)];
        t_q=(index(Number_index+1)-index(Number_index))*Ts;
        t_h_1=(index(Number_index+1)-index(Number_index-1))*Ts;
        t_h_2=(index(Number_index+2)-index(Number_index))*Ts;
        t_f_3=(index(Number_index+3)-index(Number_index-1))*Ts;
        t_f_4=(index(Number_index+4)-index(Number_index))*Ts;
        f=[f 1/5*(1/t_q/4+1/t_h_1/2+1/t_h_2/2+1/t_f_3+1/t_f_4)];
    elseif Number_index==3
        t=[t;index(Number_index) index(Number_index+1)];
        t_q=(index(Number_index+1)-index(Number_index))*Ts;
        t_h_1=(index(Number_index+1)-index(Number_index-1))*Ts;
        t_h_2=(index(Number_index+2)-index(Number_index))*Ts;
        t_f_2=(index(Number_index+2)-index(Number_index-2))*Ts;
        t_f_3=(index(Number_index+3)-index(Number_index-1))*Ts;
        t_f_4=(index(Number_index+4)-index(Number_index))*Ts;
        f=[f 1/6*(1/t_q/4+1/t_h_1/2+1/t_h_2/2+1/t_f_2+1/t_f_3+1/t_f_4)];
    elseif Number_index==L_index-1
        t=[t;index(Number_index) index(Number_index+1)];
        t_q=(index(Number_index+1)-index(Number_index))*Ts;
        t_h_1=(index(Number_index+1)-index(Number_index-1))*Ts;
        t_f_1=(index(Number_index+1)-index(Number_index-3))*Ts;
        f=[f 1/3*(1/t_q/4+1/t_h_1/2+1/t_f_1)];
    elseif Number_index==L_index-2
        t=[t;index(Number_index) index(Number_index+1)];
        t_q=(index(Number_index+1)-index(Number_index))*Ts;
        t_h_1=(index(Number_index+1)-index(Number_index-1))*Ts;
        t_h_2=(index(Number_index+2)-index(Number_index))*Ts;
        t_f_1=(index(Number_index+1)-index(Number_index-3))*Ts;
        t_f_2=(index(Number_index+2)-index(Number_index-2))*Ts;
        f=[f 1/5*(1/t_q/4+1/t_h_1/2+1/t_h_2/2+1/t_f_1+1/t_f_2)];
    elseif Number_index==L_index-3
        t=[t;index(Number_index) index(Number_index+1)];
        t_q=(index(Number_index+1)-index(Number_index))*Ts;
        t_h_1=(index(Number_index+1)-index(Number_index-1))*Ts;
        t_h_2=(index(Number_index+2)-index(Number_index))*Ts;
        t_f_1=(index(Number_index+1)-index(Number_index-3))*Ts;
        t_f_2=(index(Number_index+2)-index(Number_index-2))*Ts;
        t_f_3=(index(Number_index+3)-index(Number_index-1))*Ts;
        f=[f 1/6*(1/t_q/4+1/t_h_1/2+1/t_h_2/2+1/t_f_1+1/t_f_2+1/t_f_3)];
    else
        t=[t;index(Number_index) index(Number_index+1)];
        t_q=(index(Number_index+1)-index(Number_index))*Ts;
        t_h_1=(index(Number_index+1)-index(Number_index-1))*Ts;
        t_h_2=(index(Number_index+2)-index(Number_index))*Ts;
        t_f_1=(index(Number_index+1)-index(Number_index-3))*Ts;
        t_f_2=(index(Number_index+2)-index(Number_index-2))*Ts;
        t_f_3=(index(Number_index+3)-index(Number_index-1))*Ts;
        t_f_4=(index(Number_index+4)-index(Number_index))*Ts;
        f=[f 1/7*(1/t_q/4+1/t_h_1/2+1/t_h_2/2+1/t_f_1+1/t_f_2+1/t_f_3+1/t_f_4)];
    end
end
t=t.*1/fs;
f_GZC=sum(f)/length(f);
end