function output = bicoherence(x)
% This is file bicoherence.m. The function calculates
% the bicoherence of signals x = SP - PV
%
% References :  Choudhury et al., Automatic detection
%               and quantification of stiction in control 
%               valves (in press) + Diagnosis of poor
%               control-loop performance using
%               higher-loop statistics (brief paper)
%
%               Press et al., Numerical Recipes in Fotran 77
%
%               bicoerenza.m, matlabfile made by M.
%               Rossi
%                       
%   H. Manum (2006)

if length(x)/2 == round(length(x)/2)
else
    x = x(1:end-1);
end

% Segmentation of the input signal

LX = 64;        % segment length
LS = 128;       % DFT lengt
OV = 0.5;       % 50% overlap

k = 0;
cont = 0;
while cont == 0
	indices{k+1} = ((LX/2)*k+1:1:(LX/2)*k+LX+1);
	values = indices{k+1};
	k = k + 1;
	if values(end)>length(x) - LX/2
		cont = 1;
	end
end
numseg = length(indices);

% the while-loop above made the variable "indices" be of the correct dimention
% to use the overlapping DFT-calculation of the input signal.

% Construction of the bispectrum

for ix = 1:length(indices)
	clear mu x1 w x2 Y Bsp de1 de2
	% Remove the mean
	mu = mean(x(indices{ix}));
	x1 = x(indices{ix}) - mu;

	% Multiplication with the Hanning-window. Reference : page 547 in
	% Numerical Recipes
	for k = 1:length(x1)
		w(k) = 0.5 - 0.5*cos(2*pi*(k-1)/(LX-1));
	end

	x2 = w(:).*x1(:);

	% Calculation of the DFT (Y)
	Y = fft(x2,LS);
	for k = 1:LX
		for l = 1:LX
			Bsp(k,l) = Y(k)*Y(l)*conj(Y(k+l));
			de1(k,l) = abs(Y(k)*Y(l))^2;
			de2(k,l) = abs(Y(k+l))^2;
		end
	end
	D1{ix} = de1;	% Denominators in eq. 2 in Diagnosis... (Chudhury)
	D2{ix} = de2;
	Bs{ix} = Bsp;
end	

% Calculation of the mean of the bispectrum and the bicoherence

D1m = D1{1}; D2m = D2{1}; Bm = Bs{1};
for k = 2:length(Bs);
	D1m = D1m + D1{k};
	D2m = D2m + D2{k};
	Bm = Bm + Bs{k};
end
D1 = D1m/k; D2 = D2m/k; Bm = Bm/k;

bic2 = zeros(LX*2);
for k = 1:LX
	for l = 1:LX
		bic2(k+LX,l+LX) = abs(Bm(k,l))^2/(D1(k,l)*D2(k,l));
	end
end

for k = 1:LX*2
	for l = 1:LX*2
		if (k <= LX)&(l <= LX)
			bican(k,l) = bic2(2*LX-k+1,2*LX-l+1);
		elseif (k > LX)&(l <= LX)
			bican(k,l) = bic2(k,2*LX-l+1);
		elseif (k <= LX)&(l > LX)
		       bican(k,l) = bic2(2*LX-k+1,l);
		elseif (k > LX)&l > bic2(k,l);
            bican(k,l) = bic2(k,l);
        end 
	end
end
bicanpl = bican;
bican = bican(LX+1:LX+LS/2,LX+1:LX+LS/2);

waxis = linspace(0,0.5,length(bican)+1);
waxis = waxis(1:end-1);
limit = waxis(end)*(2/3);

k = 1;
while k < length(waxis)
	if waxis(k)>limit
		limit = k-1;
		k = Inf;
	else 
		k = k + 1;
	end
end

for f1 = 1:limit
	for f2 = 1:length(waxis)
		if f2 == 1 | f2 >= f1;
			bican(f1,f2) = 0;
		end
	end
end

for f1 = limit:length(waxis)
	for f2 = 1:length(waxis)
		if f2 == 1 | f2 >= -2*f1 + length(waxis)*2
			bican(f1,f2) = 0;
		end
	end
end

cont = 1;
for f1 = 1:length(waxis)
	for f2 = 1:min([f1, -2*f1+length(Y)]);
		bicpermed(cont) = bican(f1,f2);
		cont = cont + 1;
	end
end

bic2m = mean(bicpermed);
bic2v = std(bicpermed);
bic2max = max(bicpermed);

surf(waxis,waxis,bican');
axis([0 0.5 0 0.5 0 1]);
    %set(hcc,'view',[145 15],'Alim',[0 1],...
    %   'Clim',[0 0.08]);
    view([145 15])

[c,rows] = max(bican');
[c,column] = max(max(bican'));
f1 = rows(column);
f2 = column;
    
output.f1 = waxis(f1);
output.f2 = waxis(f2);
return

