1D cubic B-spline interpolation

Hi,

Following this discussion
http://groups-beta.google.com/group/comp.dsp/browse_frm/thread/ffc12437efe8=
4f92/73394dc71a49b15b?lnk=3Dst&q=3D&rnum=3D1&hl=3Den#73394dc71a49b15b
about which kind of interpolation suits my needs, I chose to try and
implement one-dimensional cubic B-spline interpolation.

However I failed, the result is a long way off what's supposed to
happen, and I'm afraid I understood something wrong about cubic
B-spline interpolation. Here's how I understood it : between two
original samples, you calculate new values by performing a cubic
B=E9zier interpolation using the two original samples before where your
new point is and the two other after, and in order to perform the cubic
B=E9zier interpolation you perform a series of linear interpolations as
depicted there http://en.wikipedia.org/wiki/Image:Bezier_cubic.png

Here's how I tried to implement it. The t, pX, qX and rX variables are
the same as those indicated in the link posted above. s2[i] is the
equivalent of B.

double *cubicBspline_interpolation(double *s1, int32_t n1, int32_t
factor)
{
	int32_t i, n2;
	double *s2, p0, p1, p2, p3, q0, q1, q2, r0, r1, t;

	n2 =3D n1 * factor;			//calculation of the length of the interpolated
signal
	s2 =3D malloc (n2 * sizeof(double));	//allocation of the interpolated
signal

	for (i=3Dfactor; i<n2-(factor*2); i++)	//loop going from after the two
first original samples to before the two last, i'll take care of the
particular cases once the general case will work correctly
	{
		t =3D fmod((float) i/factor, 1);	//calculation of our position between
p1 and p2
		p0 =3D s1[(i/factor)-1];		//the four p values are all samples from the
original signal
		p1 =3D s1[i/factor];
		p2 =3D s1[(i/factor)+1];
		p3 =3D s1[(i/factor)+2];
		q0 =3D p0 * (1-t) + p1 * t;	//the three q values are linearly
interpolated from the p values
		q1 =3D p1 * (1-t) + p2 * t;
		q2 =3D p2 * (1-t) + p3 * t;
		r0 =3D q0 * (1-t) + q1 * t;	//the two r values are linearly
interpolated from the q values
		r1 =3D q1 * (1-t) + q2 * t;
		s2[i] =3D r0 * (1-t) + r1 * t;	//the final value is linearly
interpolated from the r values
	}

	return s2;
}

The result you get from that is very wrong, if we can put it this way.
Filling the s2 array with either q value indicates that the linear
interpolation is correctly performed so far, however doing the same
with either r value gives a result that just doesn't look right, but I
guess it's not supposed to at this point.

I looked into documentation to make sure I understood the problem
correctly but couldn't figure out what's wrong, so I'd appreciate if
someone could explain me what I didn't get right.

Thanks in advance.