CORDUROY_ANALYSIS_util_lib.ipf

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

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Function CORD_RESCALE_execute(ListOfWaves,Source)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	string ListOfWaves
	variable Source
	
// have user enter the sample rate
		variable SampleRate = 10000, SampleFlag = 1, Xyes=1, Yyes=1, Yfactor=1, Yscale=10
		prompt Yyes, "Rescale Y axis?", popup "No;Yes"
		prompt Yfactor, "Select the scale factor (default: 10)", popup "2.5e-8;0.001;0.1;10;50;100"
		prompt Xyes, "Rescale X axis?", popup "No;Yes"
		prompt SampleFlag, "Select the sample rate (default: 10 kHz)", popup "10;25;5"
		DoPrompt "Dimension wave with correct scaling", Yyes, Yfactor, Xyes, SampleFlag
			if( V_Flag )
				return 0 								
			endif
// get the sample rate from the flag
		if(SampleFlag==2)
			SampleRate = 25000
		elseif(SampleFlag==3)
			SampleRate = 5000
		endif
// get the Y scaling factor from the flag
		switch(Yfactor)
			case 1:
				Yscale = 0.000000025
				break
			case 2:
				Yscale = 0.001
				break
			case 3:
				Yscale = 0.1
				break
			case 4:
				Yscale = 10
				break
			case 5:
				Yscale = 50
				break
			case 6:
				Yscale = 100
				break
			endswitch

// Get the delta x in ms
		variable delta = 1 / (SampleRate / 1000)
		string units = "ms"
// Set loop variables
		string CurrentWave = ""
		variable index=0
// Analysis loop starts here
		do
// Get waves from the list generated by CreateList function
			CurrentWave = StringFromList(index, ListOfWaves)
			if(!Source)
				WAVE local = TraceNameToWaveRef("",CurrentWave)			
			else
				WAVE local = $CurrentWave	
			endif
// When SelectedWaves list runs out break the loop
			if (strlen(CurrentWave) == 0)
				break
			endif
// Set the new scale for the selected wave
		switch(Xyes)
			case 2:
				SetScale/P x, 0,  1 / (SampleRate / 1000), units, local
				if(Yyes==2)
					local = local / Yscale
				endif
				break
			case 1:
				if(Yyes==2)
					local = local / Yscale
				endif
				break
		endswitch	
// Step the do loop
			index+=1
		while(1)
//		print "Rescaled Yaxis of "+ListOfWaves
end // RESCALE


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Function CORD_BASESUB_execute(ListOfWaves,Source)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	string ListOfWaves
	variable Source
	
	if(numtype(pcsr(A))==2 || numtype(pcsr(B))==2)
		variable/g abortflag=1
		Abort "The cursors must be on a trace in the top graph..."
	endif

	Variable NumWaves = ItemsInList(ListOfWaves)
	variable index=0, offset=0
	string CurrentWave = "", NewName=""
	
	string/G prefix="SUB"
	
	do

		CurrentWave = StringFromList(index, ListOfWaves)
	// When SelectedWaves list runs out break the loop
		if (strlen(CurrentWave) == 0)
			break
		endif
		
		if(!Source)
			WAVE local = TraceNameToWaveRef("",CurrentWave)			
		else
			WAVE local = $CurrentWave	
		endif

		NewName = prefix+"_"+num2str(index)
		Duplicate/O local, $NewName
		Wave local2 = $NewName
		
		offset = mean(local2,xcsr(A),xcsr(B))
		local2 -= offset

		index += 1
	while (1)

end //BaselineSubtract


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Function CORD_AVERAGE_execute(ListOfWaves,Source)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	string ListOfWaves
	variable Source
	variable NumWaves = ItemsInList(ListOfWaves), index
	string CurrentWave=""
	string DataFolder

	string/G suffix="AVG"

	do	
		CurrentWave = StringFromList(index, ListOfWaves)
		if(!Source)
			WAVE local = TraceNameToWaveRef("",CurrentWave)			
		else
			WAVE local = $CurrentWave	
		endif

	// When SelectedWaves list runs out break the loop
		if (strlen(CurrentWave) == 0)
			break
		endif

	// Generate holder average wave
		if(index==0)
			Duplicate/O local, w1
		// Get the channel name
			string ChName =  StringFromList(0, CurrentWave, "0")
		endif

	// Accumulate waves
		if(index>0)
			w1 += local
		endif
		
		index+=1
		
	while(1)
	
	w1 /= NumWaves
	string s3 = ChName+"_"+suffix
	Duplicate/O w1, $s3
	WAVE local2 = $s3
	Display/K=1 local2 as s3
	KillWaves w1
	
end // AVERAGE


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Function CORD_STDEV_execute(ListOfWaves,Source)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	string ListOfWaves
	variable Source
	variable NumWaves = ItemsInList(ListOfWaves), index
	string CurrentWave=""
	string DataFolder


	string/G suffix="SD"

	do
		
		CurrentWave = StringFromList(index, ListOfWaves)
		if(!Source)
			WAVE local = TraceNameToWaveRef("",CurrentWave)			
		else
			WAVE local = $CurrentWave	
		endif

	// When SelectedWaves list runs out break the loop
		if (strlen(CurrentWave) == 0)
			break
		elseif(index==0)
			// Get the channel name
			string ChName =  StringFromList(0, CurrentWave, "0")
			WAVE localAVG = $(":Analyzed:"+ChName+"_AVG")
				if(!WaveExists(localAVG))
					Abort "You must have an average wave first"
				endif
		endif

	// Generate holder average wave
		if(index==0)
			Duplicate/O local, w1
			w1 = 0
		endif

	// Accumulate error
		w1 += (local - localAVG)*(local - localAVG)

		index+=1
		
	while(1)

	w1 = w1 / (NumWaves-1)
	string s3 = ChName+"_"+suffix
	Duplicate/O w1, $s3
	KillWaves w1
	
end // STDEV


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Function CORD_CONCATENATE_execute(ListOfWaves,Source)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	string ListOfWaves
	variable Source
	
	string/G suffix="CAT"
	
	string dest=StringFromList(0,ListOfWaves)+"_"+suffix
	
	Concatenate/O ListOfWaves, wdest
	
	MoveWave wdest, $dest

end


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Function CORD_FILTER_execute(ListOfWaves,Source)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	string ListOfWaves
	variable Source
	
	Variable NumWaves = ItemsInList(ListOfWaves)
	variable index=0, offset=0
	string CurrentWave = "", NewName=""
	
	string/G suffix="FLT"
	string ex_str=""
	variable  filter_type=0, corner_f=0, pole_width=0
	
	prompt filter_type, "Select the filter typel",popup,"Low Pass; High Pass; Band Pass; Notch Filter"
	prompt corner_f, "Enter the corner frequency"
	prompt pole_width, "Enter either the number of poles or the appropriate width of the band pass"
	DoPrompt "Get filter information", filter_type, corner_f, pole_width
	if (V_Flag)
		return -1
	endif
		
	do

		CurrentWave = StringFromList(index, ListOfWaves)
		
		if (strlen(CurrentWave) == 0)
			break
		endif
		
		if(!Source)
			WAVE local = TraceNameToWaveRef("",CurrentWave)			
		else
			WAVE local = $CurrentWave	
		endif

		NewName = CurrentWave+"_"+suffix
		Duplicate/O local, $NewName
		Wave local2 = $NewName
		
		switch(filter_type)
			case 1:
				ex_str = "LowPass(\""+NewName+"\","+num2str(corner_f)+","+num2str(pole_width)+")"
				Execute ex_str
				break
			case 2:
				ex_str = "HighPass(\""+NewName+"\","+num2str(corner_f)+","+num2str(pole_width)+")"
				Execute ex_str
				break
			case 3:
				ex_str = "BandPass(\""+NewName+"\","+num2str(corner_f)+","+num2str(pole_width)+")"
				Execute ex_str
				break
			case 4:
				ex_str = "NotchFilter(\""+NewName+"\","+num2str(corner_f)+","+num2str(pole_width)+")"
				Execute ex_str
				break
		endswitch
		
		index += 1
		
	while (1)

End // FILTER


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Function CORD_SMOOTH_execute(ListOfWaves,Source)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	string ListOfWaves
	variable Source
	
// Prompt for the smoothing function to use
	variable SmoothType = 1, Repetitions = 1, Width=11
	prompt SmoothType, "Smoothing algorithm", popup "Binomial;Boxcar;Savitzky-Golay"
	prompt Repetitions, "Number of reps"
	prompt Width, "Box wdith"
	DoPrompt "Set smoothing parameters", SmoothType, Repetitions, Width
		if( V_Flag )
			return 0 								
		endif

	Variable NumWaves = ItemsInList(ListOfWaves)
	variable index=0
	string CurrentWave = "", NewName=""
	
	string/G suffix="SMTH"
	
// Cycle through waves calling Smooth function
	do

		CurrentWave = StringFromList(index, ListOfWaves)
		if(!Source)
			WAVE local = TraceNameToWaveRef("",CurrentWave)			
		else
			WAVE local = $CurrentWave	
		endif
		
	// When SelectedWaves list runs out break the loop
		if (strlen(CurrentWave) == 0)
			break
		endif
		
		NewName = CurrentWave+"_"+suffix
		Duplicate/O local, $NewName
		Wave local2 = $NewName	
			
		switch(SmoothType)
			case 1:
				Smooth Repetitions, local2
				break
			case 2:
				Smooth /B=(Repetitions) Width, local2
				break
			case 3:
				Smooth /S=2 Width, local2
				break
		endswitch
		
		index += 1
		
	while (1)	
	
End


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Function CORD_DIFFERENTIATE_execute(ListOfWaves,Source)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	string ListOfWaves
	variable Source

	Variable NumWaves = ItemsInList(ListOfWaves)
	variable index=0
	string CurrentWave = "", NewName=""
	
	string/G suffix = "DIF"
	
// Cycle through waves calling Differentiate function with default behaviour
	do

		CurrentWave = StringFromList(index, ListOfWaves)
		if(!Source)
			WAVE local = TraceNameToWaveRef("",CurrentWave)			
		else
			WAVE local = $CurrentWave	
		endif
		
	// When SelectedWaves list runs out break the loop
		if (strlen(CurrentWave) == 0)
			break
		endif
		
		NewName = CurrentWave+"_"+suffix
		Duplicate/O local, $NewName
		Wave local2 = $NewName
		Differentiate local2
		
		index += 1
		
	while (1)	
	
End


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Function CORD_INTEGRATE_execute(ListOfWaves,Source)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	string ListOfWaves
	variable Source
	
	Variable NumWaves = ItemsInList(ListOfWaves)
	variable index=0
	string CurrentWave = "", NewName=""
	
	string/G suffix = "INT"
	
// Cycle through waves calling Integrate function with default behaviour
	do

		CurrentWave = StringFromList(index, ListOfWaves)
		if(!Source)
			WAVE local = TraceNameToWaveRef("",CurrentWave)			
		else
			WAVE local = $CurrentWave	
		endif
		
	// When SelectedWaves list runs out break the loop
		if (strlen(CurrentWave) == 0)
			break
		endif
		
		NewName = CurrentWave+"_"+suffix
		Duplicate/O local, $NewName
		Wave local2 = $NewName
		Integrate local2
		
		index += 1
		
	while (1)
	
End

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Function CORD_CLEAN_execute(ListOfWaves,Source)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	string ListOfWaves
	variable Source
	
	Variable NumWaves = ItemsInList(ListOfWaves)
	variable index=0, jindex=0
	string CurrentWave = "", NewName=""
	
	
	
	variable threshold
	variable width=25 // in samples
	
	string/G suffix="CLN"
	
// Cycle through waves calling Integrate function with default behaviour
	do

		CurrentWave = StringFromList(index, ListOfWaves)

	// When SelectedWaves list runs out break the loop
		if (strlen(CurrentWave) == 0)
			break
		endif
		
		if(!Source)
			WAVE local = TraceNameToWaveRef("",CurrentWave)			
		else
			WAVE local = $CurrentWave	
		endif
		
		NewName = CurrentWave+"_"+suffix
		Duplicate/O local, $NewName
		Wave local2 = $NewName
		
		Differentiate local /D=tmp
		if(index==0)
			WaveStats/Q tmp
			threshold = -1*V_sdev
			print threshold
		endif
		
		// Go through and threshold the derivative to detect stim artifacts
		jindex=0
		do
			
			if(tmp[jindex] < threshold) // Use a negative threshold to avoid detecting/removing spikes
				local2[jindex-5,jindex+width] = NaN
				jindex = jindex+width
			else
				jindex+=1
			endif
			
		while(jindex<numpnts(tmp))
		
		KillWaves tmp
		
		// Linearly interpolate data to compensate for the missing data
		// Not sure how to implement this at the moment

		index += 1
		
	while (1)
	
End


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Function CORD_NEURON_Print_execute(ListOfWaves,Source)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	string ListOfWaves
	variable Source
	variable index=0
	string CurrentWave= ""

	do	
		CurrentWave = StringFromList(index, ListOfWaves)
		if (strlen(CurrentWave) == 0)	// When SelectedWaves list runs out break the loop
			break
		endif
		
		if(!Source)
			WAVE local = TraceNameToWaveRef("",CurrentWave)			
		else
			WAVE local = $CurrentWave	
		endif

		variable refNum = index+2
		Open refNum as "filename"	// Open a text file to write in data

		fprintf refNum,"%s\r",CurrentWave // First line write wavename
		fprintf refNum,"%s\r",num2str(numpnts(local))	// Second line write wave length

		variable ii					// Loop through the lines of the wave writing a new line at each point
		for(ii=0;ii<=numpnts(local);ii+=1)
			fprintf refNum,"%s\r",num2str(local[ii])
		endfor	

  	 	fprintf refNum,"%s\r"," " 	// Add a space at end of wave before beginning next wave
		index += 1					// Cycle through the waves by stepping the index
	while (1)

end


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Function CORD_MATLAB_Print_execute(ListOfWaves,Source)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// use LOAD filename -ASCII in matlab
	string ListOfWaves
	variable Source
	variable index=0
	string CurrentWave= ""

	do	
		CurrentWave = StringFromList(index, ListOfWaves)
		if (strlen(CurrentWave) == 0)	// When SelectedWaves list runs out break the loop
			break
		endif
		
		if(!Source)
			WAVE local = TraceNameToWaveRef("",CurrentWave)			
		else
			WAVE local = $CurrentWave	
		endif

		variable refNum = index+2
		Open refNum as "filename"	// Open a text file to write in data

		variable ii					// Loop through the lines of the wave writing a new line at each point
		for(ii=0;ii<=numpnts(local);ii+=1)
			fprintf refNum,"%s\r",num2str(local[ii])
		endfor	

  	 	fprintf refNum,"%s\r"," " 	// Add a space at end of wave before beginning next wave
		index += 1					// Cycle through the waves by stepping the index
	while (1)
End


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Function CORD_ASCII_Print_execute(ListOfWaves,Source)
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	string ListOfWaves
	variable Source
	variable index=0
	string CurrentWave= ""

	do	
		CurrentWave = StringFromList(index, ListOfWaves)
		if(!Source)
			WAVE local = TraceNameToWaveRef("",CurrentWave)			
		else
			WAVE local = $CurrentWave	
		endif

		if (strlen(CurrentWave) == 0)	// When SelectedWaves list runs out break the loop
			break
		endif

		variable refNum = index+2
		Open refNum as "filename"	// Open a text file to write in data

		fprintf refNum,"%s\r",CurrentWave // First line write wavename
//		fprintf refNum,"%s\r",num2str(numpnts(local))	// Second line write wave length

		variable ii					// Loop through the lines of the wave writing a new line at each point
		for(ii=0;ii<=numpnts(local);ii+=1)
			fprintf refNum,"%s\r",num2str(local[ii])
		endfor	

  	 	fprintf refNum,"%s\r"," " 	// Add a space at end of wave before beginning next wave
		index += 1					// Cycle through the waves by stepping the index
	while (1)
End


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////// UTILITY FUNCTIONS FOR THE FILTERING /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// 8/13/96 Mike Johnson (mike-johnson@llnl.gov)
// LowPass(), HighPass(), BandPass(), and NotchFilter() are digital
//   analogs of simple electronic signal filters.
// Spectrum() plots log-log frequency spectrum (amplitude or power)
//   of input data wave.
//
// LowPass() is analog of one or more identical RC low-pass filters in series
//   where corner frequency is 1/(2 Pi RC)
// Usage:  LowPass(input wave name, corner frequency, number of filters)
// Notes:  (1) Input wave assumed to have x-axis scale of time in seconds
//         (2) Output wave replaces input wave and has same number of points
//             and x-axis and precision (FFT done double precision).
//         (3) If the input wave has an odd number of points, the last point
//             is not filtered and the last 2 points of the output wave are
//             identical.
//         (4) To see a plot of filter transmission vs frequency, make a wave
//             (such as "transLowPass") whose x-axis covers the frequency
//             range of interest, then assign the y values by saying
//             transLowPass = magsqr(FLowPass(x, fZero)). See FLowPass and
//             similar complex transmission functions below in this listing.
Macro LowPass(inWave, fZero,numPole)
	Variable fZero, numPole=1
	String inWave
	
	Silent 1; PauseUpdate
	numPole = round(limit(numpole,1,Inf))
	Variable npnt = numpnts($inWave),x0 = leftx($inwave), dx = deltax($inWave)
	String unitStr = WaveUnits($inWave,0)
	If(Mod(npnt,2) != 0)	// FFT requires even number of pts
		Redimension/N=(npnt-1) $inWave
	endIf
	Variable wasSingle = 0
	If(FSinglePrecQ(inWave))  // Convert temporarily to dble precision for FFT
		Redimension/D $inWave
		wasSingle = 1
	endif
	FFT $inWave
	Variable j = 0
	Do
		$inWave *= FLowPass(x, fZero)
		j += 1
	While(j < numPole)
	IFFT $inWave; Redimension/R $inWave
	If(wasSingle)
		Redimension/S $inWave
	endif
	If(Mod(npnt,2) != 0)
		Redimension/N=(npnt) $inWave
		$inWave[npnt-1] = $inwave[npnt-2]
	endIf
	Setscale/P x,x0,dx,unitStr,$inWave
endMacro


// HighPass() is analog of one or more identical high-pass RC filters in series
//   where corner frequency is 1/(2 Pi RC)
// Usage:  HighPass(input wave name, corner frequency, number of filters)
// Notes:  (see LowPass notes)
Macro HighPass(inWave, fZero,numPole)
	Variable fZero, numPole=1
	String inWave
	
	Silent 1; PauseUpdate
	numPole = round(limit(numpole,1,Inf))
	Variable npnt = numpnts($inWave),x0 = leftx($inwave), dx = deltax($inWave)
	String unitStr = WaveUnits($inWave,0)
	If(Mod(npnt,2) != 0)	// FFT requires even number of pts
		Redimension/N=(npnt-1) $inWave
	endIf
	Variable wasSingle = 0
	If(FSinglePrecQ(inWave)) // Convert temporarily to dble precision for FFT
		Redimension/D $inWave
		wasSingle = 1
	endif
	FFT $inWave
	Variable j = 0
	Do
		$inWave *= FHighPass(x, fZero)
		j += 1
	While(j < numPole)
	IFFT $inWave; Redimension/R $inWave
	If(wasSingle)
		Redimension/S $inWave
	endif
	If(Mod(npnt,2) != 0)
		Redimension/N=(npnt) $inWave
		$inWave[npnt-1] = $inwave[npnt-2]
	endIf
	Setscale/P x,x0,dx,unitStr,$inWave
endMacro

// BandPass() is the analog of a single-resonance (Lorentzian) bandpass filter
//   (unity transmission at center frequency)
//   where the full frequency width is measured at the 1/2 power points
//   (this width definition is exact only for width << center frequency)
// Usage:  BandPass(input wave name, center frequency, full frequency width)
// Notes:  (see LowPass notes)
Macro BandPass(inWave, fZero, fWidth)
	Variable fZero, fWidth
	String inWave
	
	Silent 1; PauseUpdate
	Variable npnt = numpnts($inWave),x0 = leftx($inwave), dx = deltax($inWave)
	String unitStr = WaveUnits($inWave,0)
	If(Mod(npnt,2) != 0)	// FFT requires even number of pts
		Redimension/N=(npnt-1) $inWave
	endIf
	Variable wasSingle = 0
	If(FSinglePrecQ(inWave)) // Convert temporarily to double precision for FFT
		Redimension/D $inWave
		wasSingle = 1
	endif
	FFT $inWave
	$inWave *= FBandPass(x, fZero, fWidth)
	IFFT $inWave; Redimension/R $inWave
	If(wasSingle)
		Redimension/S $inWave
	endif
	If(Mod(npnt,2) != 0)
		Redimension/N=(npnt) $inWave
		$inWave[npnt-1] = $inwave[npnt-2]
	endIf
	Setscale/P x,x0,dx,unitStr,$inWave
endMacro

// NotchFilter() is the analog of a notch rejection filter (zero transmission
//   at center frequency, unity transmission far from center). This is the 
//   filter used in the Stanford Research SR510 lock-in amplifier for  
//   60-Hz and 120-Hz rejection.
//   The full frequency width is measured at the 1/2 power points.
//   (This width definition is exact only for width << center frequency)
// Usage:  NotchFilter(input wave name, center frequency, full frequency width)
// Notes:  (see LowPass notes)
Macro NotchFilter(inWave,fZero,fWidth)
	Variable fZero, fWidth
	String inWave
	
	Silent 1; PauseUpdate
	Variable npnt = numpnts($inWave),x0 = leftx($inwave), dx = deltax($inWave)
	String unitStr = WaveUnits($inWave,0)
	If(Mod(npnt,2) != 0)	// FFT requires even number of pts
		Redimension/N=(npnt-1) $inWave
	endIf
	Variable wasSingle = 0
	If(FSinglePrecQ(inWave)) // Convert temporarily to dble precision for FFT
		Redimension/D $inWave
		wasSingle = 1
	endif
	FFT $inWave
	$inWave *= FNotchFilter(x, fZero, fWidth)
	IFFT $inWave; Redimension/R $inWave
	If(wasSingle)
		Redimension/S $inWave
	endif
	If(Mod(npnt,2) != 0)
		Redimension/N=(npnt) $inWave
		$inWave[npnt-1] = $inwave[npnt-2]
	endIf
	Setscale/P x,x0,dx,unitStr,$inWave
endMacro

// Compute and display log-log plot of spectrum of a wave
// Two choices:  amplitude spectrum or power spectral density
// "Amplitude spectrum" means output is proportional to magnitude of
//   Fourier transform of input. Amplitude spectrum is normalized so
//   that input sine wave gives output value equal to amplitude (not 
//   peak-to-peak) of input sine wave.
// "Power spectrum" means output is proportional to magnitude squared of
//   Fourier transform of input. Power spectrum is normalized so that
//   output is spectral power density, that is, integral over all frequencies
//   of output gives mean square value of input wave.
// Usage:  Spectrum(input wave name, 1 => amplitude or 2 => power})
// Notes:  (1) Output wave name will be input wave name with suffix of 
//             "_asd" for amplitude spectrum and "_psd" for power spectrum.
//         (2) Input wave with odd number of points will have last point
//             duplicated before spectrum is taken
//         (3) Length of output wave is n/2+1, where n is length of input
//         (4) Spectrum has same precision as input wave (FFT done in 
//             double precision)
//         (5) My convention for DC value is NOT to divide by 4 in power 
//             spectrum, thus trapezoidal integration (as used by IGOR
//             "area" function) of power spectrum with non-zero DC value gives
//             correct answer, which is mean-square of input. Rectangular
//             integration (such as IGOR "Integrate") gives wrong answer for
//             non-zero DC in power spectrum.
//             My convention for DC value of amplitude spectrum is different;
//             I DO divide DC value by 2 so that DC value is mean of input.
//             (This is a big and messy subject, choose your own convention
//             if you do not like this one.)
Macro Spectrum(inWave, ampFlag)
	String inWave
	Variable ampFlag = 2
	Prompt inWave, "Compute spectrum of:"
	Prompt ampFlag, "Amplitude or power?", popup, "Amplitude;Power"
	
	Silent 1; PauseUpdate
	String outWave, Ylabel
	If(ampFlag == 1)
		outWave = inWave + "_asd"
		Ylabel = "Spectral amplitude"
	else
		outWave = inWave + "_psd"
		Ylabel = "Spectral power density [Hz\\SÐ1\\M]"
	endIf
	Variable npnt = numpnts($inWave)
	Duplicate/O $inWave, $outWave
	Variable wasSingle = 0
//	If(FSinglePrecQ(inWave)) // temporarily double precision for FFT
//		Redimension/D $outWave
//		wasSingle = 1
//	endif
	If (mod(npnt,2) != 0)
		Redimension/N=(npnt+1) $outWave
		$outWave[npnt] = $outWave[npnt-1] 
		npnt += 1
	endIf
	fft $outWave
	$outWave = cmplx(magsqr($outWave),0)  // because $outWave is complex 
	Redimension/R $outWave  // now can make $outWave real
	If(wasSingle)
		Redimension/S $outWave
	endif
	$outWave *= 2/deltax($outWave)/npnt^2  // do spectrum normalization
	If(ampFlag == 1)
		$outWave = sqrt(2*deltax($outWave)*$outWave)
		$outWave[0] /= 2
	endIf
	Display $outWave; Label left Ylabel; Label bottom "Frequency"
	Modifygraph grid=1, log = 1
endMacro	

Function/C FLowPass(x,fZero)
	Variable fZero, x
	
	Variable denom
	denom = 1+ (x/fZero)^2
	return cmplx(1/denom, x/fZero/denom)
end

Function/C FHighPass(x,fZero)
	Variable fZero, x
	
	Variable denom
	denom = 1+ (x/fZero)^2
	return cmplx((x/fZero)^2/denom, -x/fZero/denom)
end

Function/C FBandPass(x,fZero, fWidth)
	Variable fZero, x, fWidth
	
	Variable denom
	denom = (x*fWidth)^2 + (x^2-fZero^2)^2
	return cmplx((x*fWidth)^2/denom, (x^2-fZero^2)*fWidth*x/denom)
end

Function/C FNotchFilter(x,fZero, fWidth)
	Variable fZero, x, fWidth
	
	Variable factor
	factor = (x^2 - fZero^2)/((x^2 - fZero^2)^2 +(x*fWidth)^2)
	return cmplx((x^2 - fZero^2)*factor, -x*fWidth*factor)
end


////////// OLD AVERAGE FUNCTION ///////////////
//// Returns true if wave is single precision (real or complex)
//Function FSinglePrecQ(inWave)
//	String inWave
//	return (WaveType($inWave) %& 0x02)   // 0x02 means hex 02
//end
//
//// Returns true if wave is double precision (real or complex)
//Function FDoublePrecQ(inWave)
//	String inWave
//	return (WaveType($inWave) %& 0x04)   // 0x04 means hex 04
//end
//
//
//	string ListOfWaves
//	variable Source
//	variable NumWaves = ItemsInList(ListOfWaves), index
//	string CurrentWave=""
//	string DataFolder
//
//	do	
//		CurrentWave = StringFromList(index, ListOfWaves)
//		if(!Source)
//			WAVE local = TraceNameToWaveRef("",CurrentWave)			
//		else
//			WAVE local = $CurrentWave	
//		endif
//
//	// When SelectedWaves list runs out break the loop
//		if (strlen(CurrentWave) == 0)
//			break
//		elseif(index==0)
//			DataFolder = GetWavesDataFolder(local, 1)
//			SetDataFolder DataFolder
//		endif
//
//	// Make a matrix of appropriate size to hold all waves
//		if(index==0)
//			variable WaveLength = numpnts(local)
//			Make/O/N=(WaveLength, NumWaves) w1 = 0
//			Make/O/N=(NumWaves) w2 = 1
//		endif
//
//	// Assign wave to matrix column
//		w1[][index] = local
//
//		index+=1
//	while(1)
//	
//// Create standard output
//	MatrixOp/O Out0 =  w2 x w1
//	Out0 = Out0/NumWaves
//
//end // AVERAGE