2
$\begingroup$

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\)]\))"}]

enter image description here

enter image description here

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

$\endgroup$
4
  • 1
    $\begingroup$ Have you seen this? wolfram.com/xid/0pndor4xnyh79to3-bffc8e $\endgroup$
    – flinty
    Jul 7 '20 at 22:02
  • $\begingroup$ That's for PSD (for one dataset), I am looking for cross-spectrum (for two). en.wikipedia.org/wiki/Cross-spectrum $\endgroup$ Jul 7 '20 at 22:31
  • $\begingroup$ PowerSpectrumDensity can take vector time series though. It says Cross power spectral density between components... $\endgroup$
    – flinty
    Jul 7 '20 at 22:33
  • 1
    $\begingroup$ The data on that page is one vector. CROSS PowerSpectrumDensity needs two datasets. $\endgroup$ Jul 7 '20 at 22:44
1
$\begingroup$
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]}
  ]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.