Simple Spectral Analysis Functions

#pragma rtGlobals=1     // Use modern global access method.


Function CrossPower(srcX, srcY, bw_res)
    WAVE srcX, srcY
    Variable bw_res // Desired spectral resolution, Hz
   
    Variable retval
    retval = CrossPowerDensity(srcX, srcY, bw_res)
    if(retval < 0)
        return -1 // Error!
    endif
   
    WAVE/C G_XY
   
    Duplicate/C/O G_XY P_XY
    Redimension/R P_XY
    WAVE pxy = P_XY
    // Cross-power is modulus of G_XY
    pxy = cabs(G_XY)
   
    // Set units correctly; should be (srcUnits)^2 / Hz
    SetScale d, 0, 0, "(" + WaveUnits(srcX, 1) + ")^2/Hz", pxy
   
    KillWaves/Z G_XY
End

Function Coherence(srcX, srcY, bw_res, [noBiasError])
    WAVE srcX, srcY
    Variable bw_res // Desired spectral resolution, Hz
    Variable noBiasError // Flag; set to 1 to skip computation of bias error
   
    if(ParamIsDefault(noBiasError))
        noBiasError = 0
    endif
   
    // Compute Cross- and Autospectral density estimates from the source waves
    Variable retval
    retval =CrossPowerDensity(srcX, srcX, bw_res)
    if(retval < 0)
        return -1 // Error!
    endif
    WAVE/C G_XY
    Duplicate/O/C G_XY G_XX
   
    CrossPowerDensity(srcY, srcY, bw_res)
    WAVE/C G_XY
    Duplicate/O/C G_XY G_YY
   
    CrossPowerDensity(srcX, srcY, bw_res)
    WAVE/C G_XY
   
    Duplicate/O/C G_XY gamma_XY
   
    gamma_XY *= conj(G_XY) // |G_xy|^2 == G_xy * conj(G_xy)
    gamma_XY /= (G_XX * G_YY)
   
    Redimension/R gamma_XY
    // coherence is dimensionless
    SetScale d, 0, 0, "", gamma_XY
   
    Duplicate/O gamma_xy unc_gamma_xy, bias_gamma_xy
   
    NVAR Nchunk
    unc_gamma_xy = sqrt(2)*(1-gamma_xy) / (sqrt(gamma_xy) * sqrt(Nchunk)) // 1-std deviation, Bendat 11.47, p 275 (random error)
   
    if(!noBiasError)
        bias_gamma_xy = coherenceBiasError(gamma_xy, Nchunk)
    else
        KillWaves/Z bias_gamma_xy
    endif
   
    //Xcoher  = sqrt( magsqr(G_XY) / (G_XX * G_YY) )
   
    // Compute coherent output power spectrum, Bendat 3.71; g_y:x = gamma_xy(f) * g_yy(f)
    Duplicate/O G_YY, cohP_YX
    Redimension/R cohP_YX
   
    cohP_YX = gamma_XY * cabs(G_YY)
   
    // Set units correctly; should be (srcUnits)^2 / Hz (same as cross-power, since gamma_xy is dimensionless)
    SetScale d, 0, 0, "(" + WaveUnits(srcY, 1) + ")^2/Hz", cohP_YX
   
    KillWaves/Z G_XX, G_YY, G_XY
End

// coherenceBiasError: Given an estimate of coherence, and the number of averages used in the
// FFT computation of the Gxy values to derive it, returns the 1-sigma bias error.
// Equation IIe from Carter1973; IIa is available in commented-out portion due to compute inefficiency.
Function coherenceBiasError(gamSqrEst, Nd)
    Variable gamSqrEst // coherence estimated from Nd spectral averages
    Variable Nd // number of distinct FFT averages used to compute gamSqrEstimate
   
    // Equation IIa; exact, but sometimes leads to infinite compute time in hyperG2F1 function.
//  Return (1/Nd) + (Nd - 1)/(Nd + 1) * gamSqrEst * hyperG2F1(1,1,Nd+2,gamSqrEst) - gamSqrEst
    Return (1/Nd)*(1-gamSqrEst)^2 // Eq. IIe, approximation for large Nd
End

Function CrossPhase(srcX, srcY, bw_res)
    WAVE srcX, srcY
    Variable bw_res // Desired spectral resolution, Hz
   
    Variable retval
    retval = CrossPowerDensity(srcX, srcY, bw_res)
    if(retval < 0)
        return -1 // Error!
    endif
    WAVE/C G_XY
   
    Duplicate/O/C G_XY theta_XY
    Redimension/R theta_XY
    WAVE xphs = theta_XY
   
    xphs = atan2(imag(G_XY) , real(G_XY)) * 180/Pi
    // Units...
    SetScale d, 0, 0, "Degrees", xphs
   
    // Uncertainty estimate: need gamma_xy!
    coherence(srcX, srcY, bw_res, noBiasError=1)
    WAVE gamma_xy
   
    Duplicate/O xphs unc_theta_XY
    NVAR Nchunk
   
    // 1-std deviation formula; Bendat 11.62.
    unc_theta_XY = sqrt(1-gamma_xy)/(sqrt(gamma_xy)*sqrt(2*Nchunk)) * 180/Pi
   
    KillWaves/Z G_XY
End

// CrossPowerDensity: Computes the estimated cross-power spectral density G_XY given two source waves
// of identical length, units, and scaling. Creates (and overwrites if pre-existing) destination wave G_XY,
// which is the result of computing Bendat eq. 3.90, p. 72.
// Note that autopower spectral density estimates can be also computed with this function, letting the
// srcX and srcY waves be identical.
Function CrossPowerDensity(srcX, srcY, bw_res)
    WAVE srcX, srcY
    Variable bw_res // Desired spectral resolution, Hz
   
    if(numpnts(srcX) != numpnts(srcY) || !((deltax(srcX) - deltax(srcY)) < 1e-8))
        // Waves are of incompatible scale
        print "CrossPowerDensity is not yet equipped to handle differently-scaled waves. Aborting..."
        print numpnts(srcX), numpnts(srcY), numpnts(srcX)-numpnts(srcY)
        print deltax(srcX), deltax(srcY), deltax(srcX)-deltax(srcY)
        return -1
    endif
   
    fCrossPowerSpectralDensity(srcX, srcY, bw_res)
    WAVE/C G_XY
   
    // Correct Units on cross-spectral density result
    SetScale d, 0, 0, "(" + WaveUnits(srcX, 1) + ")^2/Hz", G_XY
   
End

// fCrossPowerSpectralDensity: Computes the complex cross power spectral density Gxy between waves
// w, w2 given a bandwidth resolution request bw_res.
// Expanded from the Wavemetrics fPowerSpectralDensity function, which computes a properly normalized
// (i.e. time-average power conserving) spectral density of a single wave (autopower spectral density).
// fCrossPowerSpectralDensity may be used to compute autopower spectra as well by specifying w = w2.

Function/S fCrossPowerSpectralDensity(w, w2, bw_res)
    Wave w, w2
    Variable bw_res // Desired bandwidth, Hz
//  Variable npsd               // SegLen - the number of input points in each segment
//  String windowName       // one of "Square;Hann;Parzen;Welch;Hamming;BlackmanHarris3;KaiserBessel"
    Variable removeDC   = 1
   
    Variable npsd = getNpsd(w, bw_res)
   
    String destw="G_XY"
    String srctmp=NameOfWave(w)+"_tmp"
    String srctmp2 =NameOfWave(w)+"_tm2"
    String winw=NameOfWave(w)+"_psdWin"
   
    Make/O/C/N=(npsd/2+1) $destw= 0
    WAVE/C psd= $destw
   
    Make/O/N=(npsd) $srctmp,$winw, $srctmp2
    WAVE tmp= $srctmp
    Wave tmp2 = $srctmp2
    WAVE win= $winw
    win= 1
   
    Variable winNorm
    Hanning win
    // winNorm=0.375        //  theoretical avg squared value
    WaveStats/Q win
    winNorm= V_rms*V_rms    // actual value is more accurate than a theoretical value.

                   
    Variable dc= 0
    Variable dc2 = 0
    if( removeDC )  // boolean
        WaveStats/Q w
        dc= V_Avg
        WaveStats/Q w2
        dc2 = V_Avg
    endif

    Variable psdFirst= 0
    Variable psdOffset= npsd/2
    Variable nsrc= numpnts(w)
   
    Variable nsegs
    for( nsegs= 0; psdFirst+npsd <= nsrc; nsegs += 1, psdFirst += psdOffset)
        Duplicate/O/R=[psdFirst,psdFirst+npsd-1] w, tmp
        Duplicate/O/R=[psdFirst,psdFirst+npsd-1] w2, tmp2
        if(mod(numpnts(tmp), 2)) // odd; add one point of zero padding to end of element
            InsertPoints (numpnts(tmp)-1), 1, tmp, tmp2
        endif
       
        tmp = (tmp-dc) * win
        tmp2 = (tmp2-dc) * win
        FFT/DEST=ctmp tmp               // result is a one-sided spectrum of complex values
        FFT/DEST=ctmp2 tmp2
        WAVE/C ctmp, ctmp2
        psd += conj(ctmp)*ctmp2 // summed Fourier power of one-sided spectrum
                                // (we're missing all negative frequency powers except for the Nyquist frequency)
    endfor
    Variable/G Nchunk = floor(nsegs + 1) / 2
//  Print "Used",Nchunk, "unique bins"
    CopyScales/P ctmp, psd

    // normalize the sum of PSDs
    Variable deltaF= deltax(psd)            // deltaF is the frequency bin width
    Variable norm= 2 / (npsd * npsd * nsegs * winNorm * deltaF)
    psd *= norm
    // Explanation of norm values:
    //  * 2                     total power = magnitude^2(-f) + magnitude^2(f), and we've only accumulated magnitude^2(f)
    //  / (npsd * npsd * nsegs) converts to average power
    //  / winNorm               compensates for power lost by windowing the time data.
    //  / deltaF                    converts power to power density (per Hertz)

    psd[0] /= 2         // we're not missing 1/2 the power of DC bin from the two-sided FFT, restore original value
    psd[npsd/2] /= 2    // there aren't two Nyquist bins, either.

    // Parseval's theorem (power in time-domain = power in frequency domain)
    // is satisfied if you compare:
    // time-domain average power ("mean squared amplitude" in Numerical Recipies) = 1/N * sum(t=0...N-1) w[t]^2
    // frequency-domain average power= deltaf * sum(f=0...npsd/2) destw[f]^2
   
    KillWaves/Z tmp, tmp2, win, ctmp, ctmp2
   
    return NameOfWave(psd)  // in the current data folder
End

// getNpsd: Given a source wave, computes the number of points to use for the
// 'npsd'  parameter in fCrossPowerSpectralDensity to obtain an effective bandwidth
// of bw_res.
// Also, marginally improved to improve resilience against unattainable bandwidth requests.
Function getNpsd(src, bw_res)
    WAVE src
    Variable bw_res

    NVAR/Z quiet

    Variable nsrc = numpnts(src)
    Variable nBins = (rightX(src) - leftX(src))*bw_res
    Variable nPsd = floor(4*(nsrc/(2*nBins)))
   
    Variable roundBins = floor(nBins)
    if(roundBins < 1 || nPsd > nsrc)
        nPsd = nsrc
        if(NVAR_Exists(quiet))
            if(quiet == 1)
                return nPsd
            endif
        endif
        Print "getNpsd: Requested resolution of", bw_res,"is not possible; using one bin..."
    endif
    return nPsd
End