If $\left( p_n \right)_{n=0}^{\infty}$ is a family of orthogonal polynoamials with respect to a measure $\mu$ on $[-1,1]$, and $\left( x_j, w_j \right)$ are the quadrature points and weights for the respective Gaussian quadrature rule, we can easily prove that $$ f_N (x) : \, = \sum\limits_{n=0}^{N-1} \sum\limits_{j=1}^{N} p_n (x_j) f(x_j) w_j p_n (x) \, ,$$ is the interpolation polynomial of degree $N-1$ for $f$ at the quadrature points $x_1, \ldots x_N$.
My question: While I could prove it, I couldn't find a reference for this proof in any textbook. Could you help with that?
This is cross-posted from this post in MSE.