Suppose I have a sinusoid $f(t) = A \cos(\omega t + \theta)$ and I want to evaluate it at Chebyshev points of the second kind ($\cos(\frac{2 \pi i}{N}), 0 \le i \le N, i \in \mathbb{Z}$), and then take the DFT of it to form the coefficients of an interpolating Chebyshev polynomial.
Is there any efficient way to do this beyond just doing the calculation at $N$ points and feeding that in to a FFT?
I basically want to efficiently calculate $A \cos(\omega \cos(\frac{2 \pi i}{N}) + \theta), 0 \le i \le N$, and then take the FFT of those values. Is there an identity for the embedded cosine calls or a way to predict what the DFT will look like without actually performing the FFT?
Specifically, Wikipedia gives this identity:
$\cos(t \sin(x)) = J_0(t) + 2 \sum_{k=1}^\infty J_{2k}(t) \cos(2kx)$
where $J_0$ and $J_{2k}$ are Bessel functions. The form above looks remarkably close to a DFT, which seems to suggest that FFT results are going to be something like the Bessel functions evaluated at $\omega$. But I'm a bit out of my depth.
So I had the right idea with my hunch on the Bessel functions. Found the source of the Wikipedia identity (Abramowitz and Stegun, 9.1.42-9.1.45), and found these sister identities:
$\cos( z \cos \theta) = J_0(z) + 2 \sum_{k=1}^{\infty}(-1)^kJ_{2k}(z) \cos(2k\theta)$
$\sin( z \cos \theta) = 2 \sum_{k=0}^{\infty}(-1)^kJ_{2k+1}(z)\cos((2k+1)\theta)$
The $\cos(2k\theta)$ terms look a lot like the $cos(n\pi)$ terms in a DFT. I don't entirely understand how to prove this result, but experimentation has verified that the Bessel functions evaluated at $\omega$ do indeed give the coefficients of the DFT (that is, the $\mathcal{F}_j$) of the sinusoid:
$\mathcal{F}\{a\cos( \omega \cos(\frac{2 \pi i}{N}))\}_j = \begin{cases} J_{0}(\omega) & \ j = 0\\ (-1)^{j/2}2a \cdot J_{j}(\omega) & \ j \ \text{is even}\\ 0 & \ j \ \text{is odd} \\ \end{cases}\\$
And likewise:
$\mathcal{F}\{a\sin( \omega \cos(\frac{2 \pi i}{N}))\}_j = \begin{cases} 0 & \ j \ \text{is even} \\ (-1)^{(j-1)/2}2a \cdot J_{j}(\omega) & \ j \ \text{is odd}\\ \end{cases}\\$
Importantly, calculating the first $N$ Bessel functions for the same abscissa $\omega$ is relatively efficient since most implementations of Bessel functions solve for them using recurrences anyway.
Handling the phase in the original sinusoid can be done using the identity $\cos(\alpha + \beta) = \cos \alpha \cos \beta - \sin \alpha \sin \beta$