Computing the Spectral Density Matrix of a Multivariate Time Series

Question

Paul Fishback 2012-10-30

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/52310-computing-the-spectral-density-matrix-of-a-multivariate-time-series

I'm fairly new to signal processing, and I have a multivariate time series whose spectral density matrix I wish to compute using only Fourier transform tools. For now, I'm just experimenting with a small number of channels and a small number of signals per channel. Later I'll worry about larger sets, more channels, windowing, averaging, etc.

Here's what I've tried so far, using only the FFT.

function [S]=SpectralDensityMatrixVia(Data)

% INPUT: Time series data stored in "Data" with each column corresponding to a particular channel

% OUTPUT: Spectral density matrix S, a 3-d array, where S(i,j,k) denotes the cross-spectral density between channels i and j at the k-th frequency. (I'll worry later about the appropriate scaling of the frequencies.)

[r,c]=size(Data);

nfft = 2^nextpow2(r);

% Computations for each pairwise channel combination.

for i=1:c

    for j=1:c
        xi=Data(:,i); %extract channel i readings
        xj=Data(:,j); %extract channel j readings
        Xi=fft(xi,nfft); %compute fft of channel i readings
        Xj=fft(xj,nfft); %compute fft of channel j readings
        Pij=Xi.*conj(Xj); %product of fft of i and conjugate fft of j
        for k=1:length(Xi) S(i,j,k)=Pij(k);
        end;
    end

end

Not only am I unsure as to whether this is correct, but I'm convinced there has to be an easier way, perhaps using other signal processing toolbox commands. (I understand how pwelch can be used for diagonal entries, which are power spectral densities, but what about the off-diagonal ones?)