Despite a great deal of work characterizing the statistical properties of radio frequency backscattered ultrasound signals, less is known about the statistical properties of demodulated intensity signals. Analysis of intensity is made more difficult by a strong nonlinearity that arises in the process of demodulation. This limits our ability to characterize the spatial resolution and noise properties of B-mode ultrasound images. In this paper, we generalize earlier results on two-point intensity covariance using a multivariate systems approach. We derive the mean and autocovariance function of the intensity signal under Gaussian assumptions on both the object scattering function and acquisition noise, and with the assumption of a locally shift-invariant pulse-echo system function. We investigate the limiting cases of point statistics and a uniform scattering field with a stationary distribution. Results from validation studies using simulation and data from a real system applied to a uniform scattering phantom are presented. In the simulation studies, we find errors less than 10% between the theoretical mean and variance, and sample estimates of these quantities. Prediction of the intensity power spectrum (PS) in the real system exhibits good qualitative agreement (errors less than 3.5 dB for frequencies between 0.1 and 10 cyc/mm, but with somewhat higher error outside this range that may be due to the use of a window in the PS estimation procedure). We also replicate the common finding that the intensity mean is equal to its standard deviation (i.e., signal-to-noise ratio = 1) for fully developed speckle. We show how the derived statistical properties can be used to characterize the quality of an ultrasound linear array for low-contrast patterns using generalized noise-equivalent quanta directly on the intensity signal.