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

Question

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

u1 and u2 datasets:
https://www.dropbox.com/s/shbvtsf92jg5fif/u1?dl=0\ https://www.dropbox.com/s/j1afbkffxkle08t/u2?dl=0

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?

Answer 1

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]}
  ]