The work presents a proof of convergence of the density of energy levels to a Gaussian distribution for a wide class of quadratic forms of Fermi operators. This general result applies also to quadratic operators with disorder, e.g., containing random coefficients. The spacing distribution of the unfolded spectrum is investigated numerically. For generic systems, the level spacings behave as the spacings in a Poisson process. Level clustering persists in the presence of disorder.