y=bsn2(x,N,M)

% B-spline network (BSN) filtering with dilation 2.
% for offline filtering in LFFC (learning feedforward control)
% by Dr Yangquan Chen. 17-04-2001.
% CSOIS/ECE Dept. of Utah State University
% Email: [e-mail]

% input 
%	signal x  :  eqi-spaced. 
%	length N=length(x);  
%	BSN support M  samples; SHOULD BE M=4*M4 +1
% output 
%	signal y

function y=bsn2(x,N,M)
y=x;	% store the original signal.						
% M<<N; M is odd (M-1)/2=M2; (M-1)/4=M4' M4=2*M2
% do the round off when M is even.
disp(['You specified the BSN support M as ',num2str(M)])
M4=floor((M-1)/4);
M=M4*4 +1;		% new length of support.
M2=M4*2;
disp(['Adjusted BSN support M as ',num2str(M)])
% check how many quarter B-splines.
N_hbs=floor(N/M4);
% first M4*3 points are left untouched.
% final M4*3+some points are also left untouched.
N_left=N-N_hbs*M4;
% start from point M2+1 to N-N_left-3*M4-1
for k=(3*M4+1):(N-N_left-4*M4-1)
% for k-th point, determine which 2 splines used.
n_k=floor(k/M4);	% Spline sequential number
n_p=k-n_k * M4;		% points used to determine the alpha1;
alpha1=n_p/M2;		% \in [0,1] weight for the 1st spline	
if (alpha1<0)
    disp('error here!')
    pause
end
alpha4=1-alpha1;	% weight for the 4-th spline
alpha2=0.5-alpha1;	
alpha3=1-alpha2;	% 0.5+alpha1
% range to do the averaging for spline-1,2,3,4
startp1=(n_k)*M4+1;	% for spline-1
startp2=(n_k-3)*M4+1;	% for spline-2
startp3=(n_k-1)*M4+1;	% for spline-3
startp4=(n_k-2)*M4+1;	% for spline-4
% can use a loop here, but no one can understand the algorithm ...
u1=0;
for j=1:2*M2+1;
	if (j<=M2+1)
		u1=u1+(j-1)/M2 * x(j-1+startp1);
	else
		u1=u1+(1- (j-M2-1)/M2) * x(j-1+startp1);
	end
end
u2=0;
for j=1:2*M2+1;
	if (j<=M2+1)
		u2=u2+(j-1)/M2 * x(j-1+startp2);
	else
		u2=u2+(1- (j-M2-1)/M2) * x(j-1+startp2);
	end
end
u3=0;
for j=1:2*M2+1;
	if (j<=M2+1)
		u3=u3+(j-1)/M2 * x(j-1+startp3);
	else
		u3=u3+(1- (j-M2-1)/M2) * x(j-1+startp3);
	end
end
u4=0;
for j=1:2*M2+1;
	if (j<=M2+1)
		u4=u4+(j-1)/M2 * x(j-1+startp4);
	else
		u4=u4+(1- (j-M2-1)/M2) * x(j-1+startp4);
	end
end
% the weighted output to y(.)
y(k)=(alpha1*u1+alpha2*u2+alpha3*u3+alpha4*u4)/M2/2.;
end 	%of  k