Cross spectral matrix and its Karhunen-Loeve transform

Hi all, I'm doing signal processing but have met a few problems. I'm trying to reproduce the work represented by the paper: "Tracking brain states under general anesthesia by using global coherence analysis", the algorithm is illustrated at last two pages, I simply summarize it below:

Data: x = [x1(t); x2(t); x3(t);...;xn(t)]; samples x channels, which means row vectors are time series, totally n rows (channels)

Step 1: fourier transform to each channel. obtain a power spectral matrix: P = [P1(f); P2(f); P3(f);...;Pn(f)]; each row Pn(f) is the fourier spectrum of its corresponding time-series xn(t)

Step 2: Cross-Spectral matrix

% cross spectral matrix
for i = 1: n  % n channels
    for j = i: n

          Pxx = P(i, :);
          Pyy = P(j, :);

          Pxy = sqrt(Pxx.*Pyy);

          for k = 1: nfft % frequency points of interests

              CX(i,j, k) = Pxy(k);
              CX(j,i, k) = CX(i,j, k);
          end
      end

end

CX is the cross-spectral matrix that has 3 dimensions, the 3rd dimension indicates the frequency point. E.g. CX(:,:,1) tells the cross-spectral matrix of all channels at frequency f(1).

Step 3: Karhunen-Loeve transform

P_matrix = P';     % Note that matrix P reads row-wise
R = cov(P_matrix); % covariance matrix. Matlab asks column-wise 
[V, D] = eig( R); basis = V';  % obtain orthogonal basis

% Karhunen-Loeve transform
P_ortho = basis * P;  % use row vector here!!

For a check, diag(cov(P_ortho')) should equal diag(D).

Step 4: Cross spectral matrix with orthogonal basis

Under the Karhunen-Loeve transform, the cross spectral matrix in the new basis:

for i = 1: n
    for j = i: n

          Pxx_ortho = P_ortho(i,:);
          Pyy_ortho = P_ortho(j,:);

          Pxy_ortho = real(sqrt(Pxx_ortho.*Pyy_ortho)); % not sure if it's right here

          for k = 1: nfft

              CY(i,j, k) = Pxy_ortho(k);
              CY(j,i, k) = CY(i,j, k);

          end
      end

end

The structure of CY is the same as CX, just the basis is different. Then the papers says CY is diagonal at each frequency, and this implies that the diagonal elements of CY contain the eigenvalues of CX, and the column of V contain the corresponding eigenvector at frequency f. Unfortunately, my program doesn't show this feature.

Could anyone please help to find any mistakes in my code.

Many thanks

Kyle

[Kyle]

The discuss about KLT was posted here:

http://www.mathworks.com/matlabcentral/answers/16737-confused-about-karhunen-loeve-transform