y=bsn1(x,N,M)

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

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

function y=bsn1(x,N,M)
y=x;				
% M<<N; M is odd (M-1)/2=M2
% do the round off when M is even.
disp(['You specified the BSN support M as ',num2str(M)])
M2=floor((M-1)/2);
M=M2*2 +1;		% new length of support.
disp(['Adjusted BSN support M as ',num2str(M)])

% check how many half B-splines.
N_hbs=floor(N/M2);

% first M2 points are left untouched.
% final M2+some points are also left untouched.

N_left=N-N_hbs*M2;

% start from point M2+1 to N-N_left-M2-1

for k=(M2+1):(N-N_left-M2-M2-1)

% for k-th point, determine which 2 splines used.

n_k=floor(k/M2);	% Spline sequential number
n_p=k-n_k * M2;	% points used to determine the alpha2;
alpha2=n_p/M2;		% \in [0,1] weight for the second spline	
if (alpha2<0)
    disp('error here!')
    pause
end

alpha1=1-alpha2;	% weight for the first spline

% range to do the averaging for spline-1

startp1=(n_k-1)*M2+1;
%endp1=(n_k+1)*M2 + 1;

startp2=n_k*M2 +1;
%endp2=(n_k+2)*M2 +1;

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

% the weighted output to y(.)

y(k)=(alpha1*u1+alpha2*u2)/(M2);

end 	% k