#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 Input_sound_file_(without_extension) ASP_01_shata_ISO
	word Output_file harmonic_peaks_valleys_subtracted.result 
	positive Analysis_width_(in_sec._for_Spectrogram) 0.0215332 
	integer Maximum_frequency_(Hz) 11025
	positive Step_size_(in_sec.e.g._every_10_ms) 0.01
	integer Frequency_step_(Hz) 20
	word Window_shape Hamming
	positive Vowel_onset_(in_sec._e.g._ASP_01_shata_ISO) 0.18433916150959281
	positive Vowel_offset_(in_sec._e.g._ASP_01_shata_ISO) 0.29983714946965306
endform
#subtracts harmonic valley amplitudes from harmonic peak amplitudes

filedelete 'output_file$'
select Sound 'input_sound_file$'
To Spectrogram... 'analysis_width' 'maximum_frequency' 'step_size' 'frequency_step' 'window_shape$'
spectrogramID = selected ("Spectrogram", -1)
fileappend 'output_file$' peak_sum(dB)	valley_sum(dB)	peak-valley 'newline$'

#Harmonic peak-picking starts here
step_num = 0
repeat
	spec_int = vowel_onset + step_size*step_num
	select 'spectrogramID'
	To Spectrum (slice)... 'spec_int'
	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... copy_of_spectrum
	spectrumID_copy = selected ("Spectrum", -1)
	#make every bin real-values
	Formula... if row=1 then sqrt(self^2+self[2,col]^2) else 0 fi
	#invert the spectrum
	Formula... 1/(1+self)

	select 'spectrumID'
	To Formant (peaks)... 1000
	formant_peakID = selected ("Formant (peaks)", -1)
	num_fmts = Get number of formants... 1
	sum_peak = 0
	for i from 1 to num_fmts
		peak_freq = Get value at time... 'i' 0.5 Hertz Linear

#After you get the peak frequency, identify the bin number that this peak belongs to.  Then get the band
# energy for the bin

		select 'spectrumID'
		binNum = Get bin from frequency... 'peak_freq'
		select 'spectrumCopyID'
		energy_dB = Get real value in bin... 'binNum'
###### for debugging
		fileappend 'output_file$' 'energy_dB:4' 'tab$'
		sum_peak = sum_peak + energy_dB
		select 'formant_peakID'
	endfor
	Remove

######	fileappend 'output_file$' 'sum_peak:4' 'tab$'
	fileappend 'output_file$' end_of_peaks 'tab$'
	
#Harmonic valley-picking starts here

	select 'spectrumID_copy'	
	To Formant (peaks)... 1000
	#formant peaks are actually harmonic valleys since the peaks were inverted
	formant_valleyID = selected ("Formant (peaks)", -1)
	num_fmts = Get number of formants... 1
	sum_valley = 0
	for i from 1 to num_fmts
		valley_freq = Get value at time... 'i' 0.5 Hertz Linear
		select 'spectrumID'
		binNum = Get bin from frequency... 'valley_freq'
		select 'spectrumCopyID'		
		energy_dB = Get real value in bin... 'binNum'
###### for debugging
		fileappend 'output_file$' 'energy_dB:4' 'tab$'
		sum_valley = sum_valley + energy_dB
		select 'formant_valleyID'
	endfor
	subtracted = sum_peak - sum_valley
	Remove
	select 'spectrumID'
	Remove
	select 'spectrumID_copy'
	Remove
	select 'spectrumCopyID'
	Remove
######	fileappend 'output_file$' 'sum_valley:4' 'tab$'
######	fileappend 'output_file$' 'subtracted' 'newline$'
	fileappend 'output_file$' end_of_valleys 'newline$'
	step_num = step_num + 1
#repeats until the analysis window size crosses over vowel offset
until (spec_int > vowel_offset - analysis_width)
select 'spectrogramID'
Remove