A polynomial \(F(z)=\sum_{k=0}^{n-1} \overline{z}^k \,\varphi_k(z)\) in \(\overline{z}\) with analytic coefficients \(\varphi_k(z)\) is called a polyanalytic function. One denotes by \(\mathcal{F}_p^n(\mathbb{C})\) the Fock space of all polyanalytic functions of order \(n\) such that \[ \|F\|^p=\int_{\mathbb{C}} |F(z)|^p\, e^{-\pi |z|^2/2}\, dz \] is finite. The main result provides an interpolation formula with convergence in \(\mathcal{F}_p^n(\mathbb{C})\) as follows. Denote \(\Lambda=\alpha(\mathbb{Z}+i\mathbb{Z})\) and \(\Lambda^0=\frac{1}{\alpha}(\mathbb{Z}+i\mathbb{Z})\). Set \[ S_{\Lambda^0}^n(z)=\Bigl(\frac{\pi^n}{n!}\Bigr)^{1/2} e^{\pi |z|^2} \Bigl(\frac{d}{d z}\Bigr)^n\Bigl[e^{-\pi|z|^2} \frac{(\sigma_{\Lambda^0}(z))^{n+1}}{n!z}\Bigr], \] where \(\sigma_{\Lambda^0}(z)\) denotes the classical Weierstrass sigma function associated with \(\Lambda^0\). It is proved that, if \(\alpha<1/(n+1)\), then every \(F\in \mathcal{F}_p^n(\mathbb{C})\) possesses an expansion \[ F(z)=\sum_{\lambda\in \Lambda} F(\lambda) e^{\pi\overline{\lambda} z-\pi|\lambda|^2/2} S_{\Lambda^0}^n(z-\lambda) \] interpolating its samples, with convergence in the Fock space norm. In establishing this fact, Gabor frames with Hermite functions are studied and extended to Banach frames in modulation spaces, and mapping properties of the polyanalytic Bargmann transform in modulation spaces are established.