A new technique is presented for computing noise in nonlinear circuits. The method is based on a formulation that uses harmonic power spectral densities (HPSDs), using which a block-structured matrix relation between the second-order statistics of noise within a circuit is derived. The HPSD formulation is used to devise a harmonic-balance-based noise algorithm that requires O(nN log N). time and O(nN) memory, where n represents circuit size and N the number of harmonics of the large-signal steady state. The method treats device noise sources with arbitrarily shaped PSDs (including thermal, shot and flicker noises), handles noise input correlations and computes correlations between different outputs. The HPSD formulation is also used to establish the non-intuitive result that bandpass filtering of cyclostationary noise can result in stationary noise. The new technique is illustrated using an example that exhibits noise folding and interaction between harmonic PSD components. The results are validated against Monte-Carlo simulations. The noise performance of a large industrial integrated RF circuit (with >300 nodes) is also analyzed in less than 2 hours using tile new method.