Just take the inner product of your data with a complex exponential at the frequency of interest. If $g$ is your data, then just substitute for $f$ the value of the frequency you want (e.g., 1, 3, 10, ...):
$$ \int_{-\infty}^{+\infty} g(t) e^{ -j2\pi f t } dt $$
Or if discrete:
$$ \sum_{n=0}^{N-1} g[n] e^{-j2\pi n \frac{f}{f_s} }, $$
where $f_s$ is the sampling frequency and $N$ is the number of samples you have.
The DFT (which is what the FFT calculates) implements the above but only over a subset of the frequencies. You could also interpolate the output of the FFT, but it's unnecessary and probably easier to just calculate it directly.