function [f, c, e] = cum3equalizer(y, L, Extended_CL)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          cum3equalizer.m
% This program is designed to implement Tugnait's optimization
% algorithm for SISO Equalization.
%
% f : Estimated channel response
% c : Equalizer
% y : Received signal
% L : Equalizer length
% Extended_CL : Extended channel length, which is the length of the estimated channel,
%               Usually much longer than the true one
%
% Reference: 
% [1] J.K. Tugnait, "Identification and Deconvolution of
%     Multichannel Linear Non-Gaussian Processes Using Higher Order
%     Statistics and Inverse Filter Criteria," IEEE Trans. on
%     Signal Processing, vol. 45, no. 3, pp. 658--672, March. 1997.
%
%
% Stationary Second order and third order spetrals are used.
%
% Designed by Binning Chen on October 25, 2000
%
%  Communications and Signal Processing Laboratory
%  ECE Department, Drexel University
%  Philadelphia, PA 19104, USA
%  http://www.ece.drexel.edu/CSPL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

y=y(:).';

T=length(y);

if nargin == 2
   MAX_CORRELATION_LENGTH=L;
else
   MAX_CORRELATION_LENGTH=floor(Extended_CL/2);
end


ITERATION_TIMES=30


%%% Initialization of the Equalizer
c=zeros(1, L);
c(1,floor(L/2))=1;

rho=1;

%%% Iteration BEGIN here
for iter=1:ITERATION_TIMES
   
   iter;
   
   e=filter(c(1,:),1,y(1,:));
   CUM2c=sum(e.^2)/T;
   CUM3c=sum(e.^3)/T;
   J32=CUM3c/(CUM2c^1.5);
   
   dCUM2c_dc=zeros(size(c));
   dCUM3c_dc=zeros(size(c));
   
   for jj=1:L
      dCUM2c_dc(1,jj)=e(jj:T)*y(1,1:(T-jj+1)).'*2/T;
      dCUM3c_dc(1,jj)=(e(jj:T).^2)*y(1,1:(T-jj+1)).'*3/T;
   end
   dJ32_dc=CUM2c^(-1.5)*dCUM3c_dc - 1.5*CUM3c*CUM2c^(-2.5)*dCUM2c_dc;
   dJ32_dc=dJ32_dc*sign(J32);
   
   c1=c+rho*dJ32_dc;
   c1=c1/sqrt(sum(sum(c1.^2)));
   
   if max(max(abs(c1))) < 0.001
      rho=rho/2;
      c1=c+rho*dJ32_dc;
   end
   
   %%% Calculate NEW Object function
   e=filter(c1(1,:),1,y(1,:));
   CUM2c=sum(e.^2)/T;
   CUM3c=sum(e.^3)/T;
   J32_new=CUM3c/(CUM2c^1.5);
   
   while(abs(J32_new) < abs(J32))
      rho=rho/2;
      c1=c+rho*dJ32_dc;
      
      %%% Calculate NEW Object function with reduced step size rho
      e=filter(c1(1,:),1,y(1,:));
      
      CUM2c=sum(e.^2)/T;
      CUM3c=sum(e.^3)/T;
      J32_new=CUM3c/(CUM2c^1.5);
   end
   
   J32=J32_new;
   c=c1;
   
   J32_his(iter,1)=J32;
   
end

f=xcorr(e,y,MAX_CORRELATION_LENGTH)/T;