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