#pragma TextEncoding = "Windows-1252"
#pragma rtGlobals=3		// Use modern global access method and strict wave access.

Function GeneratePulse()

Variable TimeSpan = 1 //seconds
Variable SampleRate = 10000 // Hz
Variable NrPoints = TimeSpan*SampleRate

Make/O /N=(NrPoints) TimeWave = P/SampleRate
Make/O /N=(NrPoints) DataWave


KillDataFolder/Z root:ResampledData
KilldataFolder/Z root:DifferentialData

NewDataFolder root:ResampledData
NewDataFolder root:DifferentialData

// Compute pulse as a convolution of a gaussian and exponential functions
// This "datawave" will be downsampled. 
Variable Point=1
For(Point=0; Point<NrPoints; Point+=1)
	
	Variable TimePoint=TimeWave[Point]
	variable kappa = 20 // sets slope of exponential decay
	Variable S = .2 //sets position of peak
	variable sigma = 0.05 // sets slope of gaussian rise

	Variable Term1 = kappa/sqrt(pi)
	Variable term2 = exp((kappa*(S-Timewave[Point])+((kappa*sigma)^2)/2))
	variable Term3 = 1-erf((S-timewave[point]+kappa*sigma^2)/(sqrt(2)*sigma))

	Variable result = term1*term2*term3	
	
	DataWave[Point]=Result
	
EndFor
SetScale/P x,0,(1/SampleRate), "seconds", DataWave



//Downsample the data across the entire frequency range [1, MaxFreq]
Make/O /N=(SampleRate) SampleRates=P+1
Make/O /N=(sampleRate) MaximaResampled

Display/K=1 /W=(275,0,700,250) /N=MaximaPlot root:MaximaResampled
ModifyGraph log(bottom)=1

Duplicate/O SampleRates root:ResampledData:SampledRates

Variable i
For(i=0; i<dimSize(SampleRates,0); i+=1)

	//Generate resampled wave
	Variable SamplingRate = SampleRates[i]
	Duplicate/O DataWave ResampledWave
	Resample /rate=(SamplingRate) ResampledWave
	String ResampledWaveName = "ResampledWave_"+num2str(SamplingRate)
	Duplicate ResampledWave root:ResampledData:$ResampledWaveName
	
	//For each downsampled wave find the difference between 2 consecutive data points.	
	Duplicate/O ResampledWave Diff_ResampledWave
	Variable k
	For (k=0; k<dimSize(ResampledWave,0)-1; k+=1)
		Diff_ResampledWave[k] = sqrt((ResampledWave[k] - ResampledWave[k+1])^2)	
	EndFor
	
	//Removes final data point from the Diff_ResampledWave
	If(i==dimsize(SampleRates,0)-1)
		deletePoints i+1,1,Diff_ResampledWave
	EndIf
	String DiffWaveName = "Diff_ResampledWave_"+num2str(SamplingRate)
	Duplicate/O Diff_ResampledWave root:DifferentialData:$DiffWaveName
	
	//Generate Wave Containing Maximum values of each difference between consecutive datapoints in the resampled waves.
	MaximaResampled[i] = WaveMax(Diff_ResampledWave)
	SetScale/P x,1,1, "", MaximaResampled	
	
	Print "I've done " + num2str(i+1) + " out of " + num2str(dimSize(SampleRates,0))
	
	String ListOfWindows= WinList("Diff_Plot",";","")
	If(ItemsInList(ListOfWindows)==0)
		Display/K=1 /W=(0,0,250,250) /N=Diff_Plot root:Diff_ResampledWave
		SetAxis /W=Diff_Plot left *,1		
		Else
		//SetAxis /W=Diff_Plot /A left
		//Doupdate
	EndIF
EndFor

//Find number of data points generated at each resampling frequency and plot result

setdatafolder root:ResampledData
Make/O /N=(SampleRate) NrDataPoints_W

For(k=0;k<SampleRate;k+=1)
	ResampledWaveName = "ResampledWave_" + num2str(k+1)
	Variable NrDataPoints = DimSize($ResampledWaveName,0)
	NrDataPoints_W[k] = NrDataPoints
EndFor
moveWave NrDataPoints_W root:	
setdatafolder root:
Display/K=1 /W=(0,350,250,600) /N=PlotDataPoints NrDataPoints_W vs SampleRates
SetAxis bottom *,20
SetAxis left *,20
ModifyGraph grid=1,nticks(bottom)=20,minor(bottom)=1,sep(bottom)=20
ModifyGraph mode=4,marker=8,rgb=(0,0,0)
Label left "Nr of Data Points in Wave"
Label bottom "Resampling Frequency (Hz)"
Doupdate
End