  In some applications, one is interested in reconstructing a function $f$ from
its Fourier series coefficients. The problem is that the Fourier series is
slowly convergent if the function is non-periodic, or is non-smooth. In this
paper, we suggest a method for deriving high order approximation to $f$ using a
Pad\'e-like method. Namely, by fitting some Fourier coefficients of the
approximant to the given Fourier coefficients of $f$. Given the Fourier series
coefficients of a function on a rectangular domain in $\mathbb{R}^d$, assuming
the function is piecewise smooth, we approximate the function by piecewise high
order spline functions. First, the singularity structure of the function is
identified. For example in the 2-D case, we find high accuracy approximation to
the curves separating between smooth segments of $f$. Secondly, simultaneously
we find the approximations of all the different segments of $f$. We start by
developing and demonstrating a high accuracy algorithm for the 1-D case, and we
use this algorithm to step up to the multidimensional case.
