Cross-spectral density estimates


Cross-power spectral density is performed by the Welch's averaged, modified periodogram method. The LTPDA method ao/cpsd estimates the cross-spectral density of time-series signals, included in the input aos following the Welch's averaged, modified periodogram method [1]. Data are windowed prior to the estimation of the spectra, by multiplying it with a spectral window object, and can be detrended by polinomial of time in order to reduce the impact of the border discontinuities. The window length is adjustable to shorter lenghts to reduce the spectral density uncertainties, and the percentage of subsequent window overlap can be adjusted as well.


b = cpsd(a1,a2,pl)

a1 and a2 are the 2 aos containing the input time series to be evaluated, b is the output object, and pl is an optional parameters list.


The parameter list pl includes the following parameters:

The length of the window is set by the value of the parameter 'Nfft', so that the window is actually rebuilt using only the key features of the window, i.e. the name and, for Kaiser windows, the PSLL.

As an alternative to setting the number of points 'Nfft' in each window, it's possible to ask for a given number of CPSD estimates by setting the 'Navs' parameter, and the algorithm takes care of calculating the correct window length, according to the amount of overlap between subsequent segments.

If the user doesn't specify the value of a given parameter, the default value is used.

The function makes CPSD estimates between the 2 input aos. The input argument list must contain 2 analysis objects, and the output will contain the CPSD estimate. If passing two identical objects ai, the output will be equivalent to the output of psd(ai).


The algorithm is based in standard MATLAB's tools, as the ones used by pwelch. However, in order to compute the standard deviation of mean for each frequency bin, the averaging of the different segments is performed using Welford's algorithm [2] which allows to compute mean and variance in one loop.


Evaluation of the CPSD of two time-series represented by: a low frequency sinewave signal superimposed to white noise, and a low frequency sinewave signal at the same frequency, phase shifted and with different amplitude, superimposed to white noise.

    nsecs = 1000;
    x = ao(plist('waveform','sine wave','f',0.1,'A',1,'nsecs',nsecs,'fs',10)) + ...
    y = ao(plist('waveform','sine wave','f',0.1,'A',2,'nsecs',nsecs,'fs',10,'phi',90)) + ...
    z = cpsd(x,y,plist('nfft',1000));


  1. P.D. Welch, The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms, IEEE Trans. on Audio and Electroacoustics, Vol. 15, No. 2 (1967), pp. 70 - 73
  2. B. P. Weldford, Note on a Method for Calculating Corrected Sums of Squares and Products, Technometrics, Vol. 4, No. 3 (1962), pp 419 - 420.

©LTP Team