You're almost there, you just need to be consistent in computing the discrete data values from your interpolating function. Suppose you have a continuous function $X(f)$. Let's say you want to get the discrete approximation over $N$ bins, and $f$ runs from $0$ to $f_{max}$. Then \begin{equation} X_k = \frac{1}{\Delta f} \int_{(k-1/2) \Delta f}^{(k+1/2) \Delta f} X(f) df \end{equation} You don't really have $X(f)$ of course. But you can approximate the integral on the RHS by using one of Matlab's numerical integration routines over your interpolating function.
How is this different from what you did? Well, near the 1-2 Hz peak, simply taking the downsampled DFT value misses the peak, as you mentioned. Taking the integral over the entire bin width does include this part of the signal, so you will get a lower but wider peak.
Finally, check your total signal energy. Parseval's theorem states that the Fourier transform preserves total signal energy. Consider signal points $x_i$ with $i=0 \ldots N-1$, and write the DFT as $X_i$. Then \begin{equation} \sum_{i=0}^{N-1} |x_i|^2 = \frac{1}{N} \sum_{i=0}^{N-1} |X_i|^2 \end{equation} So, any consistent resampling procedure has to keep the RHS of this equation constant. Simple interpolation and downsampling won't do the job for sure, but the procedure I've suggested should.