  (Abridged) Estimating the uncertainty on the matter power spectrum internally
(i.e. directly from the data) is made challenging by the simple fact that
galaxy surveys offer at most a few independent samples. In addition, surveys
have non-trivial geometries, which make the interpretation of the observations
even trickier, but the uncertainty can nevertheless be worked out within the
Gaussian approximation. With the recent realization that Gaussian treatments of
the power spectrum lead to biased error bars about the dilation of the baryonic
acoustic oscillation scale, efforts are being directed towards developing
non-Gaussian analyses, mainly from N-body simulations so far. Unfortunately,
there is currently no way to tell how the non-Gaussian features observed in the
simulations compare to those of the real Universe, and it is generally hard to
tell at what level of accuracy the N-body simulations can model complicated
non-linear effects such as mode coupling and galaxy bias. We propose in this
paper a novel method that aims at measuring non-Gaussian error bars on the
matter power spectrum directly from galaxy survey data. We utilize known
symmetries of the 4-point function, Wiener filtering and principal component
analysis to estimate the full covariance matrix from only four independent
fields with minimal prior assumptions. With the noise filtering techniques and
only four fields, we are able to recover the Fisher information obtained from a
large N=200 sample to within 20 per cent, for k < 1.0 h/Mpc. Finally, we
provide a prescription to extract a noise-filtered, non-Gaussian, covariance
matrix from a handful of fields in the presence of a survey selection function.
