What you need is a tensor product b-spline function, i.e.:
$$g(x,y) = \sum_{i=1}^{n_x} \sum_{j=1}^{n_y} a_{ij} B_i(x) B_i(y)$$
where $a_{ij}$ are constants to be determined by your data (for either interpolation or approximation). $n_x$ and $n_y$ are the dimensions of $a_{ij}$ matrix in each direction ($x$ or $y$).
You can solve this problem by either interpolation or approximation (recommended) using least squares:
$$\chi^2 \sim \sum_{\mu}(D^{\mu} - g(x^{\mu},y^{\mu}))^2 $$
for some data $D$ sampled from your function $f$. Although the tensor product b-spline has a 2-dimensional summation, it is easy to transform to 1-dimension and perform standard linear algebra solutions.
You will find all the theory to get you started here:
http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/index.html
You need a tensor product b-spline construction for this. Specific details in the link at units 8 and 9.
PS Tensor product is, in simplifying the definition for your needs, just a fancy way of saying that you construct a matrix $T_{ij}$ from two vectors $u_i$ and $v_j$ by constructing each matrix element $T_{ij}$ from the product:
$$T_{ij} = u_i v_j$$
and respecting the order of the indices $i,j$. Sometimes you see this written as $T_{ij} = u_i \otimes v_j$ or equivalently, $$\mathbf{T} = \mathbf{u} \otimes \mathbf{v}.$$