# How to estimates the cross power spectral density (CPSD) of two timeseries?

So, I have two sets of data that I would like to calculate the cross-spectrum for and see the phase difference between them.

This is my code for calculating the spectrum smoothed by Kaiser Window (order 3-5).

ClearAll[specWindow];
Options[specWindow] = {"Spec" -> None};

specWindow[timeSeries_, samplFreq_, overlap_, window_,
windowLen_, \[Alpha]_, OptionsPattern[]] :=

Block[{ts, tslen, ovrlp, windowlst, power, power0, freq, freq0},

If[overlap > 1., Print["overlap must be <1"]];
ts = If[EvenQ[Length[timeSeries]], timeSeries, Most[timeSeries]];
tslen = Length@ts;

ovrlp = Round[windowLen overlap, 1];
windowlst = N@Array[window[#, \[Alpha]] &, windowLen, {-1/2, 1/2}];

Switch[OptionValue["Spec"],
"Windowed",
power =
PeriodogramArray[ts, windowLen, ovrlp, windowlst][[;;
windowLen/2]];
freq =
Table[(k samplFreq)/
tslen, {k,
Range[0., tslen/2., (tslen/2. - 0.)/(windowLen/2. - 1.)]}];
Transpose[{freq, power}]
,

"Raw",
power0 = PeriodogramArray[ts][[;; tslen/2]];
freq0 = Table[(k samplFreq)/tslen, {k, Range[0., tslen/2. - 1]}];
Transpose[{freq0, power0}]
,

"Both",
power =
PeriodogramArray[ts, windowLen, ovrlp, windowlst][[;;
windowLen/2]];
freq =
Table[(k samplFreq)/
tslen, {k,
Range[0., tslen/2., (tslen/2. - 0.)/(windowLen/2. - 1.)]}];

power0 = PeriodogramArray[ts][[;; tslen/2]];
freq0 = Table[(k samplFreq)/tslen, {k, Range[0., tslen/2. - 1]}];
{Transpose[{freq, power}], Transpose[{freq0, power0}]}
]
]


The variance preserving spectra for the two datasets:

 windowLen = 3000;
specAll =
specWindow[u1, 1./\[Tau], 0.619, KaiserWindow, windowLen, 3,
"Spec" -> "Windowed"] /. {ff_, ss_} -> {3600 24. ff, ss ff};

ListLogLogPlot[specAll,
Joined -> True,
PlotRange -> {All, {10^-6, 0.2 10^-1}},
PlotStyle -> {Directive[Thickness[0.002], Darker@Blue], Red},
ImageSize -> 450, AspectRatio -> 0.5,
GridLines -> {{1.95, 5.8, 400.}, None},
FrameLabel -> {"f (CPD)",
"\!$$\*SubscriptBox[\(PSD$$, SubscriptBox[$$u$$, \
$$All$$]]\)(\!$$\*SuperscriptBox[\(m$$, $$2$$]\)\!$$\*SuperscriptBox[\ \(s$$, $$-2$$]\))"}]


Now, how can I calculate the cross-spectrum (coherence and phase) of these two?

• Have you seen this? wolfram.com/xid/0pndor4xnyh79to3-bffc8e Commented Jul 7, 2020 at 22:02
• That's for PSD (for one dataset), I am looking for cross-spectrum (for two). en.wikipedia.org/wiki/Cross-spectrum Commented Jul 7, 2020 at 22:31
• PowerSpectrumDensity can take vector time series though. It says Cross power spectral density between components... Commented Jul 7, 2020 at 22:33
• The data on that page is one vector. CROSS PowerSpectrumDensity needs two datasets. Commented Jul 7, 2020 at 22:44

fftW returns fourier transformation of data chunks with windowLen smoothed by window and overlap overlap and frequency
cutoff time defines the cutoff frequency of the LPF
Window is the LPF window, can be Kaiser, Hanning, ...
windowLen is the length of chunks
overlap is the overlap between chunks

ClearAll[fftW];
fftW[dat_, window_, windowLen_, overlap_, sampFreq_, cutofftime_] :=

Block[{freq, ovrlp, datP, datPD, datPL, fft},
Reap[

freq =
Range[0., (windowLen/2.)/2., (windowLen/4.)/(
windowLen/2. - 1.)] sampFreq/(windowLen/2.);
Sow[freq, "freq"];

ovrlp =
Round[windowLen (1 - overlap),
1];(*partition offset is counted from the start of each part,
so for say 60% overlap,
there should be a 40% offset from the start of each part*)

datP = Partition[dat, windowLen, ovrlp];
datPD = (# - Mean@#) & /@ datP;
datPL =
LowpassFilter[#, (2 \[Pi])/cutofftime 1/sampFreq, windowLen, window] & /@
datPD;(*cutofftime in sec*)
fft = Fourier /@ datPL;
Sow[fft[[All, ;; windowLen/2]], "fft"];
][[2]]
]

CPCPF gives the magnitude squared coherence, phase,  corresponding confidence bands at 95 % confidence level, and frequency

ClearAll[CPCBF];
CPCBF[v1_, v2_, \[Alpha]_, window_, windowLen_, overlap_, sampFreq_,
cutofftime_] :=

Block[{xx, yy, xy, freq, dof, x2, y2, \[Beta], picklst, coh2, coh2L, coh2U,
phase, gain, eHxy, sp, phaseU, phaseL},

xx = fftW[v1, window, windowLen, overlap, sampFreq, cutofftime][[2, 1]];
yy = fftW[v2, window, windowLen, overlap, sampFreq, cutofftime][[2, 1]];
xy = Mean[yy Conjugate@xx];

freq = fftW[v1, window, windowLen, overlap, sampFreq, cutofftime][[1, 1]];

dof = Length@xx - 1;
Print[dof];(*degree of freedom*)
\[Beta] =
1. - \[Alpha]^(
1/dof);(*I should pick coh2 values that are larger than \[Beta]*)

x2 = Mean[Abs[xx]^2];
y2 = Mean[Abs[yy]^2];

(*coherence*)
coh2 = Abs[xy]^2/(x2 y2);(*magnitude squared coherence*)

picklst = posfastGE[coh2, \[Beta]];

coh2U = Tanh[
0.5 Log[(1 + Sqrt[coh2])/(1 - Sqrt[coh2]) - 1/(2 dof) + 2.81/
Sqrt[2 dof]]]^2;
coh2L = Tanh[
0.5 Log[(1 + Sqrt[coh2])/(1 - Sqrt[coh2]) - 1/(2 dof) - 2.81/
Sqrt[2 dof]]]^2;

(*phase*)
phase = Arg /@ xy;
gain = Abs[xy]/x2;
eHxy = 2/(2 dof)
InverseCDF[FRatioDistribution[2, 2 dof], 1 - \[Alpha]] (1 - coh2) Total[
Re[Conjugate[yy] yy]]/Total[Re[Conjugate[xx] xx]];
sp = Re[ArcSin[eHxy/gain]];
phaseU = phase + sp;
phaseL = phase - sp;

{coh2[[picklst]], coh2U[[picklst]], coh2L[[picklst]], phase[[picklst]],
phaseU[[picklst]], phaseL[[picklst]], freq[[picklst]], \[Beta]};
{coh2, coh2U, coh2L, phase, phaseU, phaseL, freq, \[Beta]}
]