Cross power spectral density (CPSD) comparison between signals I want to compare various signals with Matlab's CPSD to identify shared frequency components. (If there's a better method, I would like to know!)

I have 3 signals, A, B, and C. B is pure noise while A and C share nearby frequencies.

I noticed that CPSD's output is scaled by the power input of the signal - so for comparison between A/B/C, they should have normalized power so the CPSD amplitudes are comparable to each other (hopefully?).

The issue is that if I normalize the power, the B noise signal is boosted above the true 4 Hz signal and it looks like there's minor 5 Hz stuff for the B noise signal as well - which are both false.

How can I compare A/B/C effectively and see the relations but avoid boosted noise amplitudes in CPSD output?

I've also included a periodogram which shows the stronger 5 Hz signal but not the 4 Hz signal, for reference. mscohere gave gibberish as well. Finally, I want to keep this like "exploratory" and pretend I don't know there's a low-level 4 Hz signal that I can filter for specifically.

Any suggestions would be greatly appreciated, as I could be going in a very wrong direction or doing something incorrectly and have no idea!

t = (0:6:1080)'; %sec, some time steps
rng(1863)
A = 2*sin( 1/5*t ) + 0.25*sin( 1/4*t) + 2*(rand(length(t),1)-0.5); %some data with noise  + 2*(np.random.uniform(size=t.size)-0.5)
B = 1*(rand(length(t),1)-0.5); %some data with noise
C = 0.75*sin( 1/4.8*t ) + 0.25*sin( 1/4.05*t) + 2*(rand(length(t),1)-0.5); %some data with noise

window = hamming(110); %some window
windowPad = padarray(window,size(t,1)-size(window,1),0,'post'); %padded to t size for periodogram
Fs = 1/(1); %1/sec, time delta in freq form

pwr_A = sqrt(1/length(A)*sum(A.^2)); %estimate power of signal
pwr_B = sqrt(1/length(B)*sum(B.^2)); %estimate power of signal
pwr_C = sqrt(1/length(C)*sum(C.^2)); %estimate power of signal

[Cxy_AvB,freqs_AvB] = cpsd(A,B,window,100,512,Fs);
Axy_AvB = angle(Cxy_AvB)*180/pi; 
Pxy_AvB = abs(Cxy_AvB);

[Cxy_AvC,freqs_AvC] = cpsd(A,C,window,100,512,Fs);
Axy_AvC = angle(Cxy_AvC)*180/pi; 
Pxy_AvC = abs(Cxy_AvC);

[Cxy_BvC,freqs_BvC] = cpsd(B,C,window,100,512,Fs);
Axy_BvC = angle(Cxy_BvC)*180/pi; 
Pxy_BvC = abs(Cxy_BvC);

figure(1);
subplot(2,1,1);
plot(1./freqs_AvB,Pxy_AvB);
hold on;
plot(1./freqs_AvC,Pxy_AvC);
plot(1./freqs_BvC,Pxy_BvC);
xlim([0 10]);  
xlabel('Periods (sec)');
ylabel('Arb. Power');
title('Cross power spectral density matlab');
legend('AvB(noise)','AvC','B(noise)vC', 'Location','northwest');

A = 1/pwr_A*A; %normalize to a power of 1
B = 1/pwr_B*B; %normalize to a power of 1
C = 1/pwr_C*C; %normalize to a power of 1

[Cxy_AvB,freqs_AvB] = cpsd(A,B,window,100,512,Fs);
Axy_AvB = angle(Cxy_AvB)*180/pi; 
Pxy_AvB = abs(Cxy_AvB);

[Cxy_AvC,freqs_AvC] = cpsd(A,C,window,100,512,Fs);
Axy_AvC = angle(Cxy_AvC)*180/pi; 
Pxy_AvC = abs(Cxy_AvC);

[Cxy_BvC,freqs_BvC] = cpsd(B,C,window,100,512,Fs);
Axy_BvC = angle(Cxy_BvC)*180/pi; 
Pxy_BvC = abs(Cxy_BvC);

subplot(2,1,2);
plot(1./freqs_AvB,Pxy_AvB);
hold on;
plot(1./freqs_AvC,Pxy_AvC);
plot(1./freqs_BvC,Pxy_BvC);
xlim([0 10]);  
xlabel('Periods (sec)');
ylabel('Arb. Power');
title('Normalized Power Cross power spectral density matlab');
legend('AvB(noise)','AvC','B(noise)vC', 'Location','northwest');

[pxx,f] = periodogram(A,windowPad,512,Fs);
pxx = abs(pxx);

figure(2)
plot(1./f,pxx)
hold on;

[pxx,f] = periodogram(B,windowPad,512,Fs);
pxx = abs(pxx);

plot(1./f,pxx)

[pxx,f] = periodogram(C,windowPad,512,Fs);
pxx = abs(pxx);

plot(1./f,pxx)
xlabel('Period (sec)');
ylabel('Normalized Power');
xlim([0 10]);  
title('Normalized DFT of A,B,C');
legend('A','B','C', 'Location','northwest');
I've ended up at a workable solution.

To deal with the noise, average the CPSD results of many realizations of the noise.


Pxy_AvB_mat = zeros(512/2+1,nRealizations); %preallocate
Pxy_BvC_mat = zeros(512/2+1,nRealizations); %preallocate
for i=1:nRealizations
    B = 1*(rand(length(t),1)-0.5); %new realization
    [Cxy_AvB,freqs_AvB] = cpsd(A,B,window,100,512,Fs);
    Pxy_AvB_mat(:,i) = abs(Cxy_AvB);
    [Cxy_BvC,freqs_BvC] = cpsd(B,C,window,100,512,Fs); 
    Pxy_BvC_mat(:,i) = abs(Cxy_BvC);
end
Pxy_AvB = mean(Pxy_AvB_mat,2);
Pxy_BvC = mean(Pxy_BvC_mat,2);


Use a spectral analysis tool like periodogram, pwelch, or roll your own with FFTs to verify if a signal is pure noise or not. Averaging multiple realizations of noise's spectral powers will yield a flat line - showing there are no major noise constituents. Same idea as the code above. A single realization can be misleading for spectral analysis or CPSD. Many-realizations-averaged-CPSD will still show some minor peaks with respect to the noise signal, but the spectral analysis in combination will help confirm if there are any actual frequency components (noise on its own is flat).