Formal development of the theory for harmonic analysis of neutron multiplying systems is carried out completely in the frequency domain. From basic probability theory, and an assumed reactor model, the problem is expressed as the Fokker-Planck equation in terms of the characteristic function, thus enabling the moments required for a statistical analysis to be obtained. Second-moment calculations include investigation into the bias in estimates of the power spectral density arising from the existence of finite record lengths. It is seen that for even very long records large biases can result, particularly at the lower frequencies. Variance analysis for estimates of the power spectral density investigates all moments up to and including the fourth for neutrons, delayed neutron precursors, and Fourier coefficients. The results show that for the most part, the variances can be described by a single parameter in which the extraneous neutron source plays a particularly important role. For reactors with large sources, the Fourier coefficients are shown to be Gaussian. For systems with small sources, variance in estimates of the power spectral density can become very large, and even the classical smoothed estimate is not consistent.