We present a quantum algorithm to evaluate matrix elements of functions of unitary operators. The method is based on calculating quadrature nodes and weights using data collected from a quantum processor. Given a unitary $U$ and quantum states $|ψ_0\rangle$, $|ψ_1\rangle$, the resulting quadrature rules form a functional that can then be used to classically approximate $\langleψ_1|f(U)|ψ_0\rangle$ for any function $f$. In particular, the algorithm calculates Szegö quadrature rules, which, when $f$ is a Laurent pol…