How can I compute a time-frequency cross-spectrum in MATLAB?

Question

MATLAB has a spectrogram function for the time-frequency analysis of a single signal. It also has a cpsd function for estimating the cross-frequency spectrum for two signals. However, cpsd averages across windows, collapsing the time axis into a single estimate.

Is there a function in this same family that does not average, but returns a time-frequency cross-spectrum instead?

Because I want to use this in a step-by-step tutorial for teaching purposes, I'd like to avoid the following two possible solutions:

• The cohgramc function from the Chronux toolbox (uses a multitaper spectral estimation method I don't want to get into)
• The wcoher function from the Wavelet toolbox (uses wavelets, idem)

Both of these methods return nice time-frequency cross-spectra, but I would like just to have a basic version with the familiar windowing parameters that spectrogram and cpsd use. Does this exist?

Answer 1

I don't know of a function that does all of what you ask, but it is easy enough to write. They key step is to create short-time integrations on which the FFTs may be performed. As you say, cpsd then averages all these STIs, but you can write your own to skip that step.

The key function to do this is $y = buffer (x,fftlen,overlap)$. You input your FFT processing parameters, similar to cpsd, and get a matrix of overlapping time series. I apply a Hanning window to each STI before performing FFT. Then you can simply apply $S = FFT(y)$ to get a matrix of overlapped spectra.

So it would look something like this:

% Assuming x is longer than fftlen:
overlap = fftlen/2; % 50% overlap
win = hanning(fftlen);
X = buffer(x,fftlen,fftlen/2,'nodelay');  % Matrix of overlapping STIs
numSTIs = size(X,2);
winX = X.*win(:,ones(1,numSTIs); % Time-domain windowed STIs
S = fft(winX,fftlen,1)/fftlen;  % Double-Sided Complex Spectrum Matrix
SdB = 20*log10(2*abs(S(1:fftlen/2+1,:)));  % Log Scale Single-Sided Real Spectrum Matrix

You can then calculate your auto-spectrum using either of these matices. Of course you would need an x2 to repeat and calculate a cross-spectrum matrix.

for i = 1:numSTIs
    Sx1x2(:,:,i) = S1(:,i)'*S2(:,i);
end

(If anyone knows of a way to do this without a for loop, please let me know!)