#Same as 13, but with different variable conventions
#dB value with an arbitrary overall level......based on Paul's formula
#Band energy summation for just one bin that contains the harmonic peak frequency
#Based on Paul and Yisheng's comments
#stops making spectrum if analysis window crosses vowel offset
#prompts for file names and parameters for spectrographic analysis
form Select files and parameters for spectrogram analysis
	word inputSoundFile_(without_extension) ASP_01_shata_ISO
	word outputFile harmonic_peaks_valleys_subtracted.result 
	positive analysisWidth_(in_sec._for_Spectrogram) 0.0215332 
	integer maxFrequency_(Hz) 11025
	positive stepSize_(in_sec.e.g._every_10_ms) 0.01
	integer frequencyStep_(Hz) 20
	word windowShape Hamming
	positive vowelOnset_(in_sec._e.g._ASP_01_shata_ISO) 0.18433916150959281
	positive vowelOffset_(in_sec._e.g._ASP_01_shata_ISO) 0.29983714946965306
endform
#subtracts harmonic valley amplitudes from harmonic peak amplitudes
filedelete 'outputFile$'
select Sound 'inputSoundFile$'
To Spectrogram... 'analysisWidth' 'maxFrequency' 'stepSize' 'frequencyStep' 'windowShape$'
spectrogramID = selected ("Spectrogram", -1)
fileappend 'outputFile$' sumOfPeaks(dB)	sumOfValleys(dB)	peaks-valleys 'newline$'
#Harmonic peak-picking starts here
stepNum = 0
repeat
	specLocation = vowelOnset + stepSize*stepNum
	select 'spectrogramID'
	To Spectrum (slice)... 'specLocation'
	spectrumID = selected ("Spectrum", -1)
	Copy... copyOfSpectrum
	spectrumCopyID = selected("Spectrum",-1)
	Formula... if row=1 then 10*log10(self^2+self[2,col]^2) else 0 fi
	select 'spectrumID'
	Copy... copy2OfSpectrum
	spectrumCopy2ID = selected ("Spectrum", -1)
	#make every bin real-values
	Formula... if row=1 then sqrt(self[1,col]^2+self[2,col]^2) else 0 fi
	#invert the spectrum
	Formula... 1/(1+self)

	select 'spectrumID'
	To Formant (peaks)... 1000
	formantPeaksID = selected ("Formant (peaks)", -1)
	numFormants = Get number of formants... 1
	sumPeaks = 0
	for i from 1 to numFormants
		peakFrequency = Get value at time... 'i' 0.5 Hertz Linear
		select 'spectrumCopyID'
		binNum = Get bin from frequency... 'peakFrequency'
		energy_dB = Get real value in bin... 'binNum'
###### for debugging
######		fileappend 'outputFile$' 'energy_dB:4' 'tab$'
		sumPeaks = sumPeaks + energy_dB
		select 'formantPeaksID'
	endfor
	Remove
	fileappend 'outputFile$' 'sumPeaks:4' 'tab$'
######	fileappend 'outputFile$' end_of_peaks 'tab$'
#Harmonic valley-picking starts here
	select 'spectrumCopy2ID'	
	To Formant (peaks)... 1000
	#formant peaks are actually harmonic valleys since the peaks were inverted
	formantValleysID = selected ("Formant (peaks)", -1)
	numFormants = Get number of formants... 1
	sumValleys = 0
	for i from 1 to numFormants
		valleyFrequency = Get value at time... 'i' 0.5 Hertz Linear
		select 'spectrumCopyID'
		binNum = Get bin from frequency... 'valleyFrequency'
		energy_dB = Get real value in bin... 'binNum'
###### for debugging
######		fileappend 'outputFile$' 'energy_dB:4' 'tab$'
		sumValleys = sumValleys + energy_dB
		select 'formantValleysID'
	endfor
	subtracted = sumPeaks - sumValleys
	Remove
	select 'spectrumID'
	Remove
	select 'spectrumCopyID'
	Remove
	select 'spectrumCopy2ID'
	Remove
	fileappend 'outputFile$' 'sumValleys:4' 'tab$'
	fileappend 'outputFile$' 'subtracted' 'newline$'
######	fileappend 'outputFile$' end_of_valleys 'newline$'
	stepNum = stepNum + 1
#repeats until the analysis window size crosses over vowel offset
until (specLocation > vowelOffset - analysisWidth)
select 'spectrogramID'
Remove