#############################################################################
#############################################################################
# Noise Vocoder
# Vocode all selected sounds in the Objects list
# Matthew Winn
# November 2017
#
#---------------------------------------------------------------------------
# Brief description: flexible vocoder that can be customized for specific:
# lower and upper corner frequencies, number of analysis and synthesis channels,
# an option for synthesis peak-picking (as in n-of-m / ACE cochlear implant processing),
# user-specified filter width to simulate spread of cochlear excitation, 
# user-specified cochlear shift to simulate incomplete CI insertion,
# support for customized channel-frequency allocation,
# and options for various carrier types,
# including noise, sinewaves, harmonic complexes, and low-noise noise.
# Batch processing supported for any (reasonable) number of Sound objects 
# currently open/selected in the Praat objects list. 
#---------------------------------------------------------------------------
# DISCLAIMER: I know this code isn't perfectly optimized. 
#   If you're seeking a blazing-fast (but not quite as customizable) vocoder,
#   I recommend AngelSim by Qian-Jie Fu.
#
#############################################################################
# INSTRUCTIONS 
# You fill out the vocoder parameters in the pop-up menu,
#    and it will ask you to select all the sound objects in the list 
#    that you want to vocode.
#    so, call up the sounds before you run the script. 
# It will batch-process all your selections;
#    note that some styles of vocoding are computationally intense,
#    so you will want to experiment with 1 or 2 short sounds
#    before really diving in. 
#############################################################################
# Set upper and lower frequencies for analysis filter
#----------------------------------------------------------------------------
# Set number of analysis & synthesis filters
#    (for a conventional vocoder, these will be the same number)
#    if you choose a lower number for the synthesis filter, the script
#    will create a peak-picking vocoder that emulates the 
#    advanced combination encoder (ACE)-style of peak-picking 
#    in successive time bins.
# NOTE: if you choose this option, the script will run pretty slow, as it's a serial process
#    rather than a psuedo-vectorized process. 
#    Processing time is generally twice that of the duration of the sound.
#    The conventional vocoder is *much* faster.
#----------------------------------------------------------------------------
# Choose a temporal envelope cutoff filter
#    to control temporal precision of the envelope
#----------------------------------------------------------------------------
# Choose a carrier type: noise, sinewave, or pulse-spreading harmonics
#   Pulsatile harmonic carrier is a wideband tone complex
#   whose partial phase relationships produce a pulsatile signal. 
#   pulse rate equals F0 of harmonic complex
#   (Hilkhuysen/Macherey-style high rate pulse-spreading harmonic complex
#    is in the works and will be available in the next major version)
#----------------------------------------------------------------------------
# Choose the shape of the synthesis filter as "peaked"
#    to model the spread of excitation of a cochlear implant electrode 
#    (controlled by the 'rolloff.per.mm' variable, 
#    which operates in terms of cochlear space)
# OR choose the "square" option, which is the conventional style,
#    where a flat-spectrum noise is divied up into frequency bands,
#    and the rolloff.per.mm is ignored
# if you select LNN (low-noise noise),
#    then the shape will be flat,
#    and the bandwidth will be just slightly narrower than 1 ERB
#----------------------------------------------------------------------------
# Choose the width of the synthesis filter in mm.octave
#    anything below 10 will be quite distorted, while
#    anything above 20 will usually sound pretty clear.
#    super-high numbers (> 90) will sound more like distorted sinewave vocoders
#   ** This option will only be implemented if you selected "peaked"
#      for your synthesis filter shape. ** 
#----------------------------------------------------------------------------
# Choose a length (in mm) of cochlear frequency shifting
#    to simulate incomplete insertion of the array
#    (higher numbers mean a greater upward shift in frequency;
#    shifts of 3mm and above generally make speech very difficult to understand.
#----------------------------------------------------------------------------
# If you're doing the ACE-style peak-picking, 
#    select length of the time bins over which peaks are picked. 
#    (the conventional vocoder isn't affected by this, since the envelope 
#    is carried for each channel at all times, whereas for the peak-picking style, 
#    the channels sitmulated depends on the exact part of the signal being analyzed.
#----------------------------------------------------------------------------
# You can also simulate the frequency allocation table used by the 
#   Nucleus family of devices produced by the Cochlear Corporation
# by selecting the 'simulate cochlear map' startup option.
# by selecting the Cochlear frequency table option, 
# your custom settings in the script will override your startup window settings. 
#  (note that this is NOT the actual method that Cochlear uses in their devices;
# I'm only trying to replicate some basic frequency filtering elements
#  using some basic Praat procedures.
#---------------------------------------------------------------------------
# A module has been added to simulate partial tripolar (PTP) stimulation
#    as found in the Advanced Bionics device. 
#    Note that NO attempt is made to replicate the actual processing for the AB device;
#    I'm just replicating the channel-frequency allocations and letting you choose
#    your desired settings for spread of excitation. 
#----------------------------------------------------------------------------
# If you want to simulate a specific custom frequency table,
#    you can choose that option at the bottom of the pop-up menu,
#    and it will use the channels defined at the bottom of this script.
#    It's currently set up to simulate a 15-channel cochlear implant 
#---------------------------------------------------------------------------
# By selecting the custom frequency table option, 
#    your custom settings in the script will override your startup window settings. 
#----------------------------------------------------------------------------
#  In this script, there is no attempt made to control for amplitude quantization
#     or amplitude compression.
#############################################################################

form Enter vocoder settings 
optionmenu filter 1
 option Channels based on these corner frequencies
 option use "AB PTP" map
 option use "Nucleus" map
 option use customized map
 
comment overall upper and lower corner frequencies
  real lowCornerFreq 100
  real highCornerFreq 8000

comment number of analysis & synthesis channels (n-of-m)
   integer numberStimulated 8
   integer numberOfChannels 8


comment synthesis filter parameters

optionmenu carrier_type 1
 option use noise carrier
 option use pulsed carrier
 option use sinewave
 option use PSHC
 option use original sound
 #option use TFS pulse carrier
 
real pulseRate 200

optionmenu shape: 1
   option peaked (controlled filter rolloff)
   option square
   option LNN

  real rolloff.per.mm 12

real envelope_cutoff_filter_(Hz) 600
boolean use_Hilbert 0

real scramble_mix 0.5

real compression_mult 1

# basal shift (in mm of cochlear space)
  real basal_shift_(mm) 0
  
  comment time bins for peak-picking option
   real window_refresh_rate 30


  
# run all sounds in a whole folder? 
  boolean run_folder 0
  
  sentence orig_directory L:\PraatScripts\Vocoder\Orig_sounds_demo
  sentence new_directory_name Voc_8_ch_shift_3
  boolean remove_suffix 0

endform
#############################################################################
### PROCESSES
#############################################################################
clearinfo

call setParameters

if run_folder == 0
pause select all sounds to be used for this operation
numberOfSelectedSounds = numberOfSelected ("Sound")

for thisSelectedSound to numberOfSelectedSounds
    sound'thisSelectedSound' = selected("Sound",thisSelectedSound)
endfor

for thisSound from 1 to numberOfSelectedSounds
	select sound'thisSound'
	name$ = selected$("Sound")
	print vocoding 'name$'... ('thisSound' of 'numberOfSelectedSounds') 'newline$'
	
	# as long as you're not doing the interleaved version, proceed normally
	   if keep_these_channels$ != "interleaved"
		stopwatch
		call vocode 'name$'
		time_elapsed = stopwatch
		print ...... took 'time_elapsed:3' seconds'newline$'
	   else
	   
	   call vocode_interleaved
	   
	endif
endfor

else

	Create Strings as file list: "fileList", "'orig_directory$'\*.wav"
	num_files = Get number of strings

	# only proceed if there are any files to work with
	print Processing 'num_files' files in the folder 'orig_directory$''newline$'

	if num_files < 1
	   exit no files in this directory!
	endif

	# make the new folder
	system mkdir 'orig_directory$'\'new_directory_name$'

	for file_index from 1 to num_files
	select Strings fileList

	filename$ = Get string: 'file_index'
	Read from file: "'orig_directory$'\'filename$'"
	name$ = selected$("Sound")

	call vocode 'name$'

	# save the sound
	select Sound 'name$''vocodedSuffix$'
	Save as WAV file... 'orig_directory$'\'new_directory_name$'\'name$''vocodedSuffix$'.wav

	# remove the new sound
	   select Sound 'name$''vocodedSuffix$'
	   Remove

	# remove the original sound
	   select Sound 'name$'
	   Remove

	endfor

endif

if verbose !=1
   #clearinfo
endif
call printCornerFreqs

if run_folder == 1
	# save info box to the same folder
	call saveInfoWindow 'orig_directory$'\'new_directory_name$' 'new_directory_name$'_vocoder_info
	select Strings fileList
	Remove
endif


print 'newline$'DONE!!!


#############################################################################
#############################################################################
# PROCEDURES 

# the main procedure
procedure vocode name$
   # Get Sound info
	select Sound 'name$'
	start = Get start time
	end = Get end time
	samplerate = Get sampling frequency
	
	if verbose == 1
		print vocoding 'name$''newline$'
	endif
	
   # Coerce number of signal channels variable to be 1
	numberOfSignalChannels = 1

   	# If you selected a number of synthesis channels that is smaller than
   	# the number of analysis channels, it does a peak-picking algorithm
   	# "ACE"-style processing: walk along each time bin, vocode, keep top N peaks. 
	#==============================================================================
	    # Vocode each channel the time bin
		if verbose == 1
			print creating carrier channels...'newline$'
		endif
		call createEnvelopes 'name$' numberOfChannels rolloff.per.mm envelope
		   # with all the channels made, do the peak-picking...
	     
	    # if you're peak-picking... 
	     if numberStimulated != numberOfChannels
		  if verbose == 1
			print peak-picking...'newline$'
		  endif
		  call peakpick_envelopes numberOfChannels numberStimulated window_refresh_rate
		  
		  call rename_peakpicked_envelopes
	    endif

		# make the carriers
		# create the wideband carrier that will be filtered into channels
		       if carrier_type$ = "noise"
		       	    do ("Create Sound from formula...", "carrier", numberOfSignalChannels, 0, end, samplerate, "randomGauss(0,0.1)")
			    			  
			  elsif carrier_type$ = "toneComplex"
			         # make tone complex here
			         Create Sound as tone complex: "carrier", 0, end, samplerate, "Sine", pulseRate, 0, 0, 0
			  
			  elsif carrier_type$ = "sinewave"
			  	# ... need to wait until channel loop
			   
			  elsif carrier_type$ = "PSHC"
			   	  call open_PSHC .duration .samplerate
			  elsif carrier_type$ == "original_sound"
			   	# this will require a new per-channel process
			   	pause original-sound carrier not ready yet, as it requires per-channel filtering first
			   	    select Sound 'name$'
			   	    # make TFS
			   	    call do_Hilbert_transform '.sound$'
			   	    select Sound 'name$'_ENV
			   	    Remove
			   	    select Sound 'name$'_TFS
			   	    Rename... carrier
			endif
		
			   
	# create carrier channels for each frequency band
	for .thisChannel from 1 to numberOfChannels
		    if verbose = 1
		    	print creating carrier for channel '.thisChannel''newline$'
		     endif
		if carrier_type$ = "sinewave"
			cf = carrier_freq_center_'.thisChannel'
		   	call make_sinewave cf end samplerate carrier
		endif
		
		# if it's LNN that you want, make it now,
		# and do the envelope multiplication inline
		if shape$ = "LNN"
		   cf = carrier_freq_center_'.thisChannel'
		   call make_LNN_ERB cf end samplerate carrier
		   	
		# if it's not LNN, multiply the carrier by the envelope now
		 endif
    
	   	    # multiply the noise by the intensity tier from the original sound channel
			select Sound carrier
			Copy... broad_carrier_channel_'.thisChannel'_w_ENV
			Formula... self[col] * Sound_channel_'.thisChannel'_ENV[col]
		
	   	 # Filter the output channel
	   	 
	   	 call filter_carrier_channel .thisChannel
	   	 
	   	 if carrier_type$ = "sinewave"
	   	 	select Sound carrier
	   	 	Remove
	   	 endif
		# end loop through channels
	endfor

	# combine the channels
	   select Sound carrier_channel_1
	   Copy... temp_vocoded
	   if numberOfChannels > 1
	   for chan_index from 2 to numberOfChannels
		Formula... self[col] + Sound_carrier_channel_'chan_index'[col]
	   endfor

	call cleanup_temp_objects

	#call re_tilt_spectrum

	call finalize_vocode_process
endproc

# # ## ### ##### ######## ############# #####################
# # ## ### ##### ######## ############# 
# # ## ### ##### ######## 
# # ## ### ##### 
# # ## ### 
# # ## 
# # 
#
# PROCEDURES 
#

procedure saveInfoWindow outputDirectory$ outputFileName$
	# save all the contents of the info window
	filedelete 'outputDirectory$'\'outputFileName$'.txt
	fappendinfo 'outputDirectory$'\'outputFileName$'.txt
endproc

procedure make_LNN_ERB cf .duration .samplerate .newname$
  num_iterations = num_LNN_iterations
  
  # make a broadband noise
  do ("Create Sound from formula...", "noise_to_filter", 1, 0, .duration, .samplerate, "randomGauss(0,0.1)")
  
  # filter it into a single ERB  
  erb_center = hertzToErb(cf)
  erb_width = erb(cf)
  if verbose == 1
  	#print 'cf:0' Hz is 'erb_center:3' ERB'newline$'
  	#print 'cf:0' Hz has ERB of 'erb_width:3''newline$'
  endif
  
  .freq_lower = cf - (erb_width/2) + 5
  .freq_upper = cf + (erb_width/2) - 5
    
  skirt = 10
  target_intensity = 60
  
  # divide by envelope, repeat filtering  
  # Make low-noise noise
  
  original_sound$ = "noise_to_filter"
  
  # initiate a starting sound
     select Sound 'original_sound$'
     Filter (pass Hann band): .freq_lower, .freq_upper, skirt
     Rename... 'original_sound$'_LNN_0
     Scale intensity... target_intensity
  
  for n from 1 to num_iterations
  if verbose == 1
  	if n < 2
  	   print LNN iteration number 'n'
  	else
  	   print ... 'n'
  	endif
  	if n == num_iterations
   	   print 'newline$'
  	endif
  endif
  m = n - 1
     # which sound to analyze & multiply?
     # use the last output of this process
     # (if it's run #1, use the original copy,
     # which was renamed to end in "_0")
     
     target_sound$ = "'original_sound$'_LNN_'m'"
  
          # get the envelope and fine structure
  	if use_Hilbert == 1
  		call get_envelope_TFS 'target_sound$' 8000
  		select Sound 'target_sound$'_TFS
  		Remove
  	else
		# maybe this is a place to sub in the alternate envelope extraction
		# so that it doesn't use Hilbert. 
		# have separate Hilbert-LNN
		# and a Praat-LNN options
		call fast_envelope 'target_sound$'
  	endif
  	
  	select Sound 'target_sound$'
  	
  	Copy... temp
  	# divide by its own envelope
  	   Formula: "self[col] / Sound_'target_sound$'_ENV[col]"
  	
  	# filter out spectral splatter
  	   Filter (pass Hann band): .freq_lower, .freq_upper, skirt
  
  	# rename for the next iteration
  	Rename... 'original_sound$'_LNN_'n'
  	
  	Scale intensity... target_intensity
  
  	select Sound temp
  	plus Sound 'target_sound$'_ENV
  	Remove
  
  endfor 
 
   # cleanup the steps along the way
	select Sound noise_to_filter
	Remove
	for n from 1 to num_iterations
	    index_to_remove = n - 1
	    select Sound 'original_sound$'_LNN_'index_to_remove'
	    Remove
	endfor

	# rename new output sound
	select Sound 'original_sound$'_LNN_'num_iterations'
	Rename... '.newname$'

endproc

procedure make_sinewave .freq .duration .samplerate .newname$
	Create Sound from formula: "'.newname$'", 1, 0, .duration, .samplerate, "1/2 * sin(2*pi*.freq*x)"
endproc

procedure get_envelope_TFS .name$ .lpf
  
    select Sound '.name$'
  
  # 1: Time-domain to frequency-domain conversion (DFT)
    spectrum = To Spectrum: "no"
    Rename: "original"
  
  # 2: Hilbert transform
    spectrumHilbert = Copy: "hilbert"
    Formula: "if row=1 then Spectrum_original[2,col] else -Spectrum_original[1,col] fi"
    soundHilbert = To Sound
  
  # 3: Obtain the ENV from the analytic signal
    env = Copy: "'.name$'_ENV"
    Formula: "sqrt(self^2 + Sound_'.name$'[]^2)"
    
    # low-pass filtered version
    if .lpf > 0
       Rename... temp
       Filter (pass Hann band)... 0 .lpf 10
       Rename... '.name$'_ENV
       select Sound temp
       Remove
    endif
  
    
  # 4: Obtain the TFS (method 1: cosine of the angle of the analytic signal)
    selectObject: soundHilbert
    tfs = Copy: "'.name$'_TFS"
    Formula: "cos(arctan2(self, Sound_'.name$'[]))"
    
    
  # 5: cleanup
    removeObject: spectrum, spectrumHilbert, soundHilbert
  
endproc

procedure modify_ENV .name$ env_LPF .method$ .keep_periodicity .scramble_mix .suffix$ .draw
  # constants
   #compression = 0.5
   env_skirt = 10
   #.draw = 0
   cleanup_comparison = 1
   .clean_mix = 1 - .scramble_mix

   #--------------------------------#
    selectObject: "Sound '.name$'"

  # get some basic info
	rms = Get root-mean-square: 0, 0
	envelope_intensity = Get intensity (dB)
	.samplerate = Get sampling frequency

  # obtain the periodicity envelope 
     # (everything above the LF envelope cutoff)
     # highpass the modified envelope to get periodicity spectrum
	Filter (pass Hann band): env_LPF, 600, env_skirt
	Rename... '.name$'_periodicity_envelope
	#print ... periodicity 

   # lowpass to get the slow envelope
	selectObject: "Sound '.name$'"
	Filter (pass Hann band): 0, env_LPF, env_skirt
	Rename: "'.name$'_envelope_LPF"

   # modify the LF envelope

   if .method$ == "scramble"
   # print using SCRAMBLE method to modify envelope'newline$'
    # randomize phase of the LPF envelope
	call randomize_spectrum_phase '.name$'_envelope_LPF env_LPF
	Rename... '.name$'_envelope_LPF_modified
	
	print ... phase proc done

   elsif .method$ == "invert"
      # print using INVERT method to modify envelope'newline$'
      # invert envelope
 	  call invert_sound '.name$'_envelope_LPF
	  print ... inverted envelope
	select Sound '.name$'_envelope_LPF_inverted
	Rename... '.name$'_envelope_LPF_modified
	
   else
       pause method$ '.method$' not recognized! 
   endif

#    # make sure envelope is all positive
# plax correct to zero here or later?? 
	call correct_to_zero '.name$'_envelope_LPF_modified

   # add the periodicity envelope
    # i.e. add low band with complementary high band 
	select Sound '.name$'_envelope_LPF_modified
	Copy... '.name$'_envelope_LPF_modified_w_periodicity
	
   if .keep_periodicity == 1
	Formula: "self[col] + Sound_'.name$'_periodicity_envelope[col]"
	print ... periodicity added
   endif

    # fix intensity of the envelope
	Scale intensity... envelope_intensity

   # add a constant so that it's all positive
	Copy... '.name$'_envelope_LPF_modified_w_periodicity_corrected
	call correct_to_zero '.name$'_envelope_LPF_modified_w_periodicity_corrected

   # name the final output
   	Rename... to_resample
   	Resample... .samplerate 20
   	
	Rename... '.name$''.suffix$'
	
   # Mix it with the original clean mix
   	Formula... self[col]*.scramble_mix + Sound_'.name$'[col]*'.clean_mix'
   		
	select Sound to_resample
	Remove

   if .draw = 1
	call draw_envelopes '.name$'_envelope_LPF '.name$'_envelope_LPF_modified
   endif

   # cleanup
	select Sound '.name$'_periodicity_envelope
	if cleanup_comparison = 1
	   plus Sound '.name$'_envelope_LPF
	   plus Sound '.name$'_envelope_LPF_modified
	endif
	plus Sound '.name$'_envelope_LPF_modified_w_periodicity
	
	if .method$=="scramble"
	   plus Sound Phases
	   plus Spectrum Ranphas
	   plus Spectrum '.name$'_envelope_LPF
	endif
   Remove
   select Sound '.name$''.suffix$'
endproc


procedure randomize_spectrum_phase .name$ .lpf

    ## the following lines are from Paul Boersma's script at http://uk.groups.yahoo.com/group/praat-users/message/59 	
	select Sound '.name$'
	orig_duration = Get total duration
	.samplerate = Get sampling frequency
	
	downsample_envelope = 1
	if downsample_envelope = 1
	# downsample
	   Resample... env_LPF*4 50
	   Rename... temp
	
	   select Sound '.name$'
	   Remove
	
	   select Sound temp
	   Rename... '.name$'
	   print ... downsampled 
	endif
	To Spectrum... yes

	select Spectrum '.name$'
	nbin = Get number of bins
	   # report number of bins
	   # (it gets really small!)
	   # pause number of bins is 'nbin'
	
	Create Sound... Phases 0 1 nbin randomUniform (0, 2 * pi)

   #You then create a new Spectrum object with the same power
       # as the original, but with a scrambled phase:

	select Spectrum '.name$'
	Copy... Ranphas
	Formula... if row = 1
	... then Spectrum_'.name$' [1, col] * cos (Sound_Phases [1, col])
	... + Spectrum_'.name$' [2, col] * sin (Sound_Phases [1, col])
	... else Spectrum_'.name$' [1, col] * sin (Sound_Phases [1, col])
	... + Spectrum_'.name$' [2, col] * cos (Sound_Phases [1, col])
	... fi
	
	print ... phase randomized 

   # sonify the spectrum
        print about to sonify... 
	To Sound
	Rename... temp
	print ... sonified
	
	# upsample if you downsampled
	   if downsample_envelope = 1
	      Resample... .samplerate 50
	      Rename... temp2
	      
	      select Sound temp
	      Remove
	      select Sound temp2
	      Rename... temp
	      print ... upsampled
	   endif

	Extract part: 0, orig_duration, "rectangular", 1, "no"
	Rename... '.name$'_scrambled_phase

	select Sound temp
	Remove
	select Sound '.name$'_scrambled_phase
endproc


procedure draw_envelopes .sound1$ .sound2$
	Erase all
    # figure out how high the y axis should go 
	select Sound '.sound2$'
	max_peak = Get maximum: 0, 0, "None"
	select Sound '.sound1$'
	orig_peak = Get maximum: 0, 0, "None"

	if orig_peak > max_peak
	   max_peak = orig_peak
	endif

	Line width: 1
	Red
	select Sound '.sound2$'
	Draw: 0, 0, 0, max_peak, "no", "Curve"
	# draw original envelope
	Line width: 2
	Black
	select Sound '.sound1$'
	Draw: 0, 0, 0, 0, "yes", "Curve"
	Line width: 1
endproc

procedure correct_to_zero .sound$
   select Sound '.sound$'
	min = Get minimum: 0, 0, "None"
	if min < 0 
	Formula: "self[col] - min"
	endif
endproc

procedure invert_sound .sound$
	select Sound '.sound$'
	Copy... '.sound$'_inverted
	.max_peak = Get maximum: 0, 0, "None"
	Formula... .max_peak - self[col]
endproc

procedure divideIntoAnalysisChannels .name$
   # this divides the *original* sound into frequency channels
   # 
   # Pre-emphasis has **already happened**
   # so that picked peaks are not biased toward low-frequency channels
   	    select Sound '.name$'
   	    Copy... '.name$'_to_analyze
	
   # extract each analysis channel using a sequence of FFT filters
      # with 50 Hz skirt
     for .thisChannel from 1 to numberOfChannels
	#print Creating channel '.thisChannel' analysis channel 'newline$'
	low = analysis_freq_low_'.thisChannel'
	high = analysis_freq_high_'.thisChannel'
	width = 50

	select Sound '.name$'_to_analyze
	Filter (pass Hann band)... low high width
	select Sound '.name$'_to_analyze_band
	Rename... '.name$'_channel_'.thisChannel'
     endfor
     
   # cleanup
        select Sound '.name$'_to_analyze
        Remove
endproc

procedure filterMMRolloff .name$ .cf .rolloff.per.mm .suffix$
   # this procedure filters the carrier channel 
   # to have a rolloff in terms of dB/mm in cochelar space 
   # Greenwood cochlea parameters set elsewhere in the script
     select Sound '.name$'
     Filter (formula)...  if x > 1 
		... then self*10^(-(abs((log10((x/aA)+k)*length/a)-(log10((cf/aA)+k)*length/a))*.rolloff.per.mm)/20) else self fi
     Rename... '.name$''.suffix$'
endproc

procedure createEnvelopes .sound$ .numAnalysisChannels .rolloff.per.mm .envelopeCutoffFreq
   # This is the workhorse function series. 
   # First, get some basic info
       select Sound '.sound$'
       .duration = Get total duration
       .samplerate = Get sampling frequency
       
    #if control.for.tilt == 1
    	# first preemphasize
    	Filter (pre-emphasis): 200
    	Rename... '.sound$'
    endif
    
    # divide analysis sound into analysis channels
    select Sound '.sound$'
    
	call divideIntoAnalysisChannels '.sound$'
	#yields: '.sound$'_channel_'.thisChannel' for all channels

    # Remove the temporarily preemphasized sound
    # because it shares a common name with the original sound object
    #if control.for.tilt == 1
    	
    	select Sound '.sound$'
    	Remove
    endif
    
	# MAKE THE ENVELOPES
    for .thisChannel from 1 to .numAnalysisChannels
	   # first get the original intensity   
		select Sound '.sound$'_channel_'.thisChannel'
		channel_'.thisChannel'_intensity = Get intensity (dB)
		if channel_'.thisChannel'_intensity < 0 
		   channel_'.thisChannel'_intensity = 0.01
		endif
	   
	   if use_Hilbert == 1
	   # get Hilbert transform envelope
	     call do_Hilbert_transform '.sound$'_channel_'.thisChannel'
	      # envelope is "'.sound$'_channel_'.thisChannel'_ENV"
	      # TFS is "'.sound$'_channel_'.thisChannel'_TFS"
	   	# remove TFS (we don't need that)
	   	select Sound '.sound$'_channel_'.thisChannel'_TFS
	   	Remove
	   else
	   	# do fast envelope using Intensity Tier, etc
	   	call fast_envelope '.sound$'_channel_'.thisChannel'
	   	# envelope is "'.sound$'_channel_'.thisChannel'_ENV"
	   	
	   
	   # low-pass filter the envelope
		if .envelopeCutoffFreq > 0 
			select Sound '.sound$'_channel_'.thisChannel'_ENV
			# rename for easy object removal later
			Rename... full_envelope
			
		     # just extract the low-pass envelope (i.e. below 50 Hz)
			Filter (pass Hann band): 0, .envelopeCutoffFreq, 2
			Rename... envelope_LPF
			if verbose == 1
				print envelope has been low-pass filtered
	   		endif
			
			select Sound full_envelope
			Remove
			
		elsif .envelopeCutoffFreq == 0
		    # scramble envelope phase
			select Sound '.sound$'_channel_'.thisChannel'_ENV
			env_LPF = 50
			keep_periodicity = 1
			draw_env = 0
			call modify_ENV '.sound$'_channel_'.thisChannel'_ENV env_LPF scramble keep_periodicity scramble_mix _scrambled draw_env
			select Sound '.sound$'_channel_'.thisChannel'_ENV_scrambled
			Rename... envelope_LPF
			
		elsif .envelopeCutoffFreq == -1
		    # invert the envelope
			select Sound '.sound$'_channel_'.thisChannel'_ENV
			env_LPF = 50
			keep_periodicity = 1
			draw_env = 0
			call modify_ENV '.sound$'_channel_'.thisChannel'_ENV env_LPF scramble keep_periodicity scramble_mix _inverted draw_env
			select Sound '.sound$'_channel_'.thisChannel'_ENV_inverted
			Rename... envelope_LPF
		endif
		
		if verbose == 1
			if .thisChannel < 2
			print envelopes have been low-pass filtered at '.envelopeCutoffFreq' Hz 'newline$'
			endif
		endif
		
		
	
	   # compress the envelope
	   	call compress_sound envelope_LPF 1-compression_mult
	   	# yields envelope_LPF_compressed
	   	
	   	if verbose == 1
	   		print envelope has been compressed
	   	endif
	   	
	   	select Sound envelope_LPF_compressed
	   	Rename... channel_'.thisChannel'_ENV
	   	
	    # cleanup
		select Sound envelope_LPF
		#plus Sound envelope_LPF_compressed
		plusObject ("Sound '.sound$'_channel_'.thisChannel'")
		if envelope <= 0 
		   plus Sound '.sound$'_channel_'.thisChannel'_ENV
		endif
		#plusObject ("Sound carrier_channel_'.thisChannel'")
		Remove
	endfor
		
endproc


procedure fast_envelope .name$
   # make the amplitude envelope
	select Sound '.name$'
	.samplerate = Get sampling frequency
	.duration = Get total duration
	
	To Intensity: 800, 0, "yes"
	selectObject: "Intensity '.name$'"
	Down to IntensityTier
	To AmplitudeTier
	Down to TableOfReal
	To Matrix
	Transpose
	To Sound (slice): 2
	Rename... temp
	Scale times to: 0, .duration
	Resample: .samplerate, 5
	Rename: "'.name$'_ENV"

   #cleanup
	select Intensity '.name$'
	plus IntensityTier '.name$'
	plus AmplitudeTier '.name$'
	plus TableOfReal '.name$'
	plus Matrix '.name$'
	plus Matrix '.name$'_transposed
	plus Sound temp
	Remove
endproc

procedure peakpick_envelopes num_channels num_selected window_refresh_rate
	# yields
	# "channel_'n'_ENV_peakpicked"

	# debug default values
	#num_channels = 8
	#num_selected = 4
	#window_refresh_rate = 30
	num_channels_dropped = num_channels - num_selected
	
	selectObject: "Sound channel_1_ENV"
	orig_samplerate = Get sampling frequency

	# grab all the channel envelopes... 
	   selectObject: "Sound channel_1_ENV"
	   for chan_index from 2 to num_channels
	      plusObject: "Sound channel_'chan_index'_ENV"
	   endfor

	# filter to anti-alias, then downsample
		Filter (pass Hann band): 0, window_refresh_rate, 20
		resamplerate = window_refresh_rate*2
		resamplerate = 'resamplerate:0'
		Resample: resamplerate, 20

	# correct to fix everything >= 0
	    for chan_index from 1 to num_channels
		#select Sound channel_'chan_index'_ENV_band_'resamplerate'
		#minimum = Get minimum: 0, 0, "None"
		#if minimum  < 0
		#Formula: "self[col] - 'minimum'"
		#endif
		# just eliminate values below 0
		Formula: "if self[col] < 0 then 0 else self[col] endif"
	    endfor


	select Sound channel_1_ENV_band_'resamplerate'
	num_samples = Get number of samples

	Create simple Matrix: "xy", num_channels, num_samples+1, "0"
	num_rows = Get number of rows
	num_cols = Get number of columns

	# channel loop
	for row_index from 1 to num_rows
	   # set channel number in first column
	   Set value: row_index, 1, row_index

	   for col_index from 2 to num_cols
	      # set channel number
	       Set value: row_index, col_index, Sound_channel_'row_index'_ENV_band_'resamplerate'[col_index-1]
	   endfor
	endfor

	select Matrix xy
	To TableOfReal
	To Table: "junk"
	select Matrix xy
	plus TableOfReal xy
	Remove

	select Table xy
	Remove column: "junk"
	Set column label (index): 1, "channel"

	for col_index from 2 to num_cols
	   col_sample_index = col_index-1
	   Set column label (index): col_index, "s'col_sample_index'"
	endfor

	# loop through each envelope sample
	for col_index from 2 to num_cols
		sample_index = col_index - 1
		select Table xy
		# sort by intensity in this time window sample
		Sort rows: "s'sample_index'"

		# loop through dropped channels
		for n_dropped from 1 to num_channels_dropped
			# zero-out the dropped channel
			Set numeric value: n_dropped, "s'sample_index'", 0
		endfor
	endfor

	# re-sort by channel
		select Table xy
		Sort rows: "channel"

	# remove channel column
		Remove column: "channel"
		Down to Matrix

	for chan_index from 1 to num_channels
		select Matrix xy
		To Sound (slice): chan_index
		Override sampling frequency: resamplerate
		Shift times to: "start time", 0

		Formula... if self[col] > 0 then 1 else 0 endif
		Rename... junk
		Resample: orig_samplerate, 2
		Rename: "channel_'chan_index'_ENV_peakpicked_downsampled_binary"
		select Sound junk
		Remove
	endfor

	# must now multiply the original envlope (which contained periodicity)
	# with the peak-picked envelope
	# so that the periodicity is maintaned,
	# but the entire envelope energy is affected
	for chan_index from 1 to num_channels
		# multiple binary peak-picked envelope
		# with true envelope
		select Sound channel_'chan_index'_ENV
		Copy... channel_'chan_index'_ENV_peakpicked
		Formula... self[col] * Sound_channel_'chan_index'_ENV_peakpicked_downsampled_binary[col]
	endfor

	# major cleanup of all things apart from peakpicked channel envelopes

	select Matrix xy
	plus Table xy
	for n from 1 to num_channels
	   plus Sound channel_'n'_ENV_band
	   plus Sound channel_'n'_ENV_band_'resamplerate'
	   plus Sound channel_'n'_ENV_peakpicked_downsampled_binary
	endfor
	Remove

	# now ready to re-multiply peakpicked envelope and TFS

endproc

procedure rename_peakpicked_envelopes
  for n from 1 to numberOfChannels
     select Sound channel_'n'_ENV_peakpicked
     select Sound channel_'n'_ENV
     Remove
     select Sound channel_'n'_ENV_peakpicked
     Rename... channel_'n'_ENV
  endfor
endproc

procedure filter_carrier_channel .thisChannel
     if shape$ =="peaked"
	    # FILTER using formula centered around [carrier_freq_center_'thisChannel' ]
	    # this is where "current spread" is simulated 
			if verbose == 1
				print filtering sound with peaked filter'newline$'
			endif
		select Sound broad_carrier_channel_'.thisChannel'_w_ENV
		cf = carrier_freq_center_'.thisChannel'
		call filterMMRolloff broad_carrier_channel_'.thisChannel'_w_ENV cf rolloff.per.mm _filt
		Rename... carrier_channel_'.thisChannel'

	elsif shape$ =="square"
		if verbose == 1
			print filtering sound with square filter'newline$'
		endif
	   # Band-pass filter using corner frequencies that define a flat-spectrum rectangular channel
		select Sound broad_carrier_channel_'.thisChannel'_w_ENV
		call bandpass broad_carrier_channel_'.thisChannel'_w_ENV analysis_freq_low_'.thisChannel' analysis_freq_high_'.thisChannel' 50 carrier_channel_'.thisChannel'

	elsif shape$=="LNN"
		# here's where we should do an LNN procedure.
		# but first filter into a single ERB 
		# centered at the center frequency of the channel. 
		freq_high = make_LNN_ERB.freq_upper
		freq_low = make_LNN_ERB.freq_lower
		call bandpass broad_carrier_channel_'.thisChannel'_w_ENV freq_low freq_high 50 carrier_channel_'.thisChannel'
		
		select Sound carrier
		Remove
	elsif carrier_type$ == "sinewave"
		# for sinwaves, process them as if they were LNNs
		cf = carrier_freq_center_'.thisChannel'
		erb_center = hertzToErb(cf)
		  erb_width = erb(cf)
		  
		  .freq_lower = cf - (erb_width/2) + 5
		  .freq_upper = cf + (erb_width/2) - 5
		call bandpass broad_carrier_channel_'.thisChannel'_w_ENV .freq_lower .freq_upper 50 carrier_channel_'.thisChannel'

	endif
	select Sound carrier_channel_'.thisChannel'
	Scale intensity... channel_'.thisChannel'_intensity
endproc


procedure evenOrOdd .which$ .number
   # .which$ is the choice of whether "odd" or "even" 
   # is the class assigned to 'evenOrOdd.keep'
	.even = '.number' mod 2
	if '.even' <>0
		## number is odd
		odd = 1
		even = 0
		if .which$ = "odd"
			   .keep = 1
			else
			   .keep = 0
		endif
	else
		odd = 0
		even = 1
		if .which$ = "even"
			   .keep = 1
			else
			   .keep = 0
		endif
	endif
	if .which$ = "all"
		.keep = 1
	endif
endproc



procedure vocode_interleaved
   # if you're doing an interleaved-channel verson, 
   # then make the odd, then even, then combine them
	keep_these_channels$ = "odd"
	call vocode 'name$'
	select Sound 'name$''vocodedSuffix$'
	Rename... 'name$''vocodedSuffix$'_left

	keep_these_channels$ = "even"
	call vocode 'name$'
	select Sound 'name$''vocodedSuffix$'
	Rename... 'name$''vocodedSuffix$'_right

	plus Sound 'name$''vocodedSuffix$'_left
	Combine to stereo
	Rename... 'name$''vocodedSuffix$'

	# intensity 
	Scale intensity... orig_intensity

	# correct for clipping
		peak = do ("Get absolute extremum...", 0, 0, "None")
		if peak > 0.99
		   do ("Scale peak...", 0.99)
		endif

	# cleanup L and R channels 
	select Sound 'name$''vocodedSuffix$'_right
	plus Sound 'name$''vocodedSuffix$'_left
	Remove
endproc

procedure setChannelCornerFrequencies lowCornerFreq highCornerFreq numberOfChannels
	## Calculate upper & lower cochlear position boundaries
	## using the inverse Greenwood function
		lowCornerPos = log10(('lowCornerFreq'/'aA')+'k')*'length'/'a'
		highCornerPos = log10(('highCornerFreq'/'aA')+'k')*'length'/'a'
	
    ### Set LOW corner cochlear positions for each channel 
    	## create the first low corner frequency
    	# directly from the frequency on the input form 
		loPos1 = lowCornerPos
	# establish the cochlear position of the LOW-freq boundary
	# of the analysis channels
		for thisChannel from 2 to (numberOfChannels+1)
		   prevChannel = thisChannel-1

	   # calculate space interval and add it to the previous landmark 
		   loPos'thisChannel' = ((highCornerPos - lowCornerPos)/(numberOfChannels))+loPos'prevChannel'
		endfor

        # Set HIGH corner cochlear POSITIONS for each channel
	for thisChannel from 1 to numberOfChannels
		nextChannel = thisChannel +1
		
		hiPos'thisChannel' = loPos'nextChannel'
	
	# Set CENTER positions 
	# (ideal electrode positions) for each channel 
		centerPos'thisChannel' = ((hiPos'thisChannel'-loPos'thisChannel')/2) + loPos'thisChannel'

	endfor

    ## Set the variable for carrier FREQUENCIES
	   for thisChannel from 1 to numberOfChannels

		analysis_freq_low_'thisChannel' 	= 'aA'*((10^('a'*loPos'thisChannel'/'length'))-'k')
		analysis_freq_center_'thisChannel' 	= 'aA'*((10^('a'*centerPos'thisChannel'/'length'))-'k')
		analysis_freq_high_'thisChannel' 	= 'aA'*((10^('a'*hiPos'thisChannel'/'length'))-'k')

	   # set electrode position (carrier) in the frequency domain, 
	   # including place shift
		carrier_freq_center_'thisChannel' = 'aA'*((10^('a'*(centerPos'thisChannel'+shift)/'length'))-'k')
	   endfor
	
   ## Now you have established lots of global variables that you can call in other parts of the script
endproc

procedure bandpass .sound$ .lf .hf .width .newname$
   # This is just a wrapper for the native Praat function
   # written as a one-liner to include object renaming. 
   select Sound '.sound$'
   Filter (pass Hann band): .lf, .hf, .width
   Rename... '.newname$'
endproc

procedure printCornerFreqs
   # print header
	print Channel'tab$'Analysis Low'tab$'Analysis Center'tab$'Analysis High'tab$'Cochlear position'tab$'Shift'tab$'Final Cochlear Position'newline$'

   # print subsequent lines 
    for thisChannel from 1 to numberOfChannels
      # establish temporary variables for printing to the info window
      		loFreqAnalysis = analysis_freq_low_'thisChannel'
		centFreqAnalysis = analysis_freq_center_'thisChannel'
		hiFreqAnalysis = analysis_freq_high_'thisChannel'
		
		centerPos = centerPos'thisChannel'
		
		final.cochlear.position = centerPos + shift
	print 'thisChannel''tab$''loFreqAnalysis:0''tab$''centFreqAnalysis:0''tab$''hiFreqAnalysis:0''tab$''centerPos:2''tab$''shift''tab$''final.cochlear.position:2''newline$'

    endfor
endproc

procedure cleanup_temp_objects
	select Sound channel_1_ENV
	if carrier_type$ != "sinewave"
		plus Sound carrier
	endif

	#plus Sound broad_carrier_channel_1
	plus Sound broad_carrier_channel_1_w_ENV
	if remove_vocoded_channel == 1
	   plus Sound carrier_channel_1
	endif

	if numberOfChannels > 1
	   for chan_index from 2 to numberOfChannels
		plus Sound channel_'chan_index'_ENV
		#plus Sound broad_carrier_channel_'chan_index'
		plus Sound broad_carrier_channel_'chan_index'_w_ENV
		if remove_vocoded_channel == 1
		   plus Sound carrier_channel_'chan_index'
		endif
	   endfor
	endif
	Remove
endproc

procedure tiltSpectrum .name$ .lowFiltRange .highFiltRange .lowFreqMult .highFreqMult .newname$
	select Sound '.name$'
	.samplerate = Get sampling frequency
	Filter (formula)... if x>.lowFiltRange and x<=.highFiltRange 
	... then self * 10^ (((((x-.lowFiltRange)/(.highFiltRange-.lowFiltRange)) * (.highFreqMult-.lowFreqMult) + .lowFreqMult)/ 20)) else self fi
	
	Filter (formula)... if x>.highFiltRange 
	... then self * 10^((((('.samplerate'/2) - x)/(('.samplerate'/2)-'.highFiltRange'))*'.highFreqMult')/20) else self fi
	Rename... '.newname$'
	
	select Sound '.name$'_filt
	Remove

	select Sound '.newname$'
endproc

procedure re_tilt_spectrum
pause what does this sound like before re-tilting?
	if control.for.tilt = 1
	   # if you did pre-emphasis, do de-emphasis now 
	   # to return to original spectral shape
  	   # and then restore the original name ("temp")
		call tiltSpectrum temp_vocoded analysis_freq_center_1 analysis_freq_center_'numberOfChannels' 0 -6 temp_vocoded_tilt
		select Sound temp_vocoded
		Remove
		select Sound temp_vocoded_tilt
	     # remove low-frequency artifact
		Filter (stop Hann band)... 0 50 40
		Rename... temp_vocoded
		select Sound temp_vocoded_tilt
		Remove
	endif
endproc

procedure finalize_vocode_process
	# add it to the zeroed original signal,
	# to ensure equal time domains
		select Sound 'name$'
		orig_intensity = Get intensity (dB)
		Copy... 'name$''vocodedSuffix$'
		Formula... Sound_temp_vocoded [col]
		
	# de-emphasize here, if you pre-emphasized
		#if control.for.tilt == 1
			Filter (de-emphasis): 200
			Rename... temp1
			Filter (stop Hann band)... 0 70 40
			Rename... temp
			select Sound 'name$''vocodedSuffix$'
			plus Sound temp1
			Remove
			select Sound temp
			Rename... 'name$''vocodedSuffix$'
		endif
		
	# ramp the onset & offset
		call ramp_onset_offset 'name$''vocodedSuffix$' 3
		# rename bc of name change
		select Sound 'name$''vocodedSuffix$'_ramped
		Copy... temp
		select Sound 'name$''vocodedSuffix$'_ramped
		plus Sound 'name$''vocodedSuffix$'
		Remove
		select Sound temp
		Rename... 'name$''vocodedSuffix$'

	# scale to original intensity
		Scale intensity... orig_intensity

	# correct for clipping
		peak = do ("Get absolute extremum...", 0, 0, "None")
		if peak > 0.99
			do ("Scale peak...", 0.99)
		endif

	#==============================================================================
	select Sound temp_vocoded
	Remove
endproc



procedure ramp_onset_offset .name$ .ramptime_ms
	# convert seconds to milliseconds
	.ramptime = .ramptime_ms/1000
	
	select Sound '.name$'
	Copy... '.name$'_ramped
	start1 = Get start time
	end1 = Get end time

	    Formula... if x<start1 + .ramptime  
	       ...then self * ((x-start1)/.ramptime) 
	       ...else self endif  

	     Formula... if x>('end1' - '.ramptime')  
		...then self * (1-((x-end1 + '.ramptime')/'.ramptime'))
		...else self endif
	endif
endproc


procedure open_PSHC .duration .samplerate
   # This wont work for anybody except me, sorry folks :)
   # it is an experimental method to use a pulse-spreading
   # harmonic complex as a carrier rather than noise/sine tones. 
   
	# first find pulse rate used in PSHC code (round down)
		pulseRate.f0 = pulseRate div 100

	# open up the PSHC rather than create it from scratch 
	#Read from file: "C:\Users\Matt\Documents\PraatScripts\Vocoder\pshc_5_1_3200_10_500.wav"
	Read from file: "C:\Users\Matt\Documents\PraatScripts\Vocoder\pshc_'pulseRate.f0'_1_12000_10_'pulseRate'.wav"
	Rename... temp
	# extract the duration meeded to match that of the vocoded sound. 
	Extract part: 0, .duration, "rectangular", 1, "no"
	Resample: .samplerate, 50
	Rename... temp_resampled
	selectObject: "Sound temp"
	plusObject: "Sound temp_part"
	Remove
	selectObject: "Sound temp_resampled"
	# rename it "noise",
	# since that's what the subsequent functions expect. 
	Rename... noise
endproc

procedure do_Hilbert_transform .name$

  select Sound '.name$'

# 1: Time-domain to frequency-domain conversion (DFT)
  spectrum = To Spectrum: "no"
  Rename: "original"

# 2: Hilbert transform
  spectrumHilbert = Copy: "hilbert"
  Formula: "if row=1 then Spectrum_original[2,col] else -Spectrum_original[1,col] fi"
  soundHilbert = To Sound

# 3: Obtain the ENV from the analytic signal
  env = Copy: "'.name$'_ENV"
  Formula: "sqrt(self^2 + Sound_'.name$'[]^2)"

  
# 4: Obtain the TFS (method 1: cosine of the angle of the analytic signal)
  selectObject: soundHilbert
  tfs = Copy: "'.name$'_TFS"
  Formula: "cos(arctan2(self, Sound_'.name$'[]))"
  
  
# 5: cleanup
  removeObject: spectrum, spectrumHilbert, soundHilbert

endproc

procedure compress_sound .sound$ .comp_ratio
	# compress the envelope
	# assumes all values are between 0 and 1
	# compressed_value = (((1 - value) * (1 - comp_ratio)) + (comp_ratio) ) * value
	# you can manually set max peak to be 1,
	# but this compresses it relative to its own peak. 
	   selectObject: "Sound '.sound$'"
	   #max_peak = Get absolute extremum: 0, 0, "None"
	   #Copy: "'.sound$'_compressed"
	   #Formula: "(((max_peak - abs(self[col])) * (max_peak - .comp_ratio)) + (.comp_ratio) ) * self[col]"

	# formula from Praat Vocal Toolkit
	
	wrk = selected("Sound")
		pow1 = Get power... 0 0
		intensity = To Intensity... 400 0 0
		Formula... -self*'.comp_ratio'
		intensitytier = Down to IntensityTier
		plus wrk
		tmp = Multiply
		pow2 = Get power... 0 0
		max = Get maximum... 0 0 None
		nocheck Scale...  max*sqrt(pow1/pow2)
		Rename... '.sound$'_compressed
		
		select intensity
		plus intensitytier
		#plus tmp
		
		Remove

endproc


procedure simulate_cochlear
   # match channel corner frequencies for the Nucleus device
   # as indicated on the clinical mapping software. 
	numberOfChannels = 22
	numberStimulated = 8
	
	analysis_freq_low_1  = 188
	analysis_freq_low_2  = 313
	analysis_freq_low_3  = 438
	analysis_freq_low_4  = 563
	analysis_freq_low_5  = 688
	analysis_freq_low_6  = 813
	analysis_freq_low_7  = 938
	analysis_freq_low_8  = 1063
	analysis_freq_low_9  = 1188
	analysis_freq_low_10 = 1313
	analysis_freq_low_11 = 1563
	analysis_freq_low_12 = 1813
	analysis_freq_low_13 = 2063
	analysis_freq_low_14 = 2313
	analysis_freq_low_15 = 2688
	analysis_freq_low_16 = 3063
	analysis_freq_low_17 = 3563
	analysis_freq_low_18 = 4063
	analysis_freq_low_19 = 4688
	analysis_freq_low_20 = 5313
	analysis_freq_low_21 = 6063
	analysis_freq_low_22 = 6938
		
	analysis_freq_center_1  = 250
	analysis_freq_center_2  = 375
	analysis_freq_center_3  = 500
	analysis_freq_center_4  = 625
	analysis_freq_center_5  = 750
	analysis_freq_center_6  = 875
	analysis_freq_center_7  = 1000
	analysis_freq_center_8  = 1125
	analysis_freq_center_9  = 1250
	analysis_freq_center_10 = 1425
	analysis_freq_center_11 = 1650
	analysis_freq_center_12 = 1925
	analysis_freq_center_13 = 2175
	analysis_freq_center_14 = 2500
	analysis_freq_center_15 = 2875
	analysis_freq_center_16 = 3300
	analysis_freq_center_17 = 3800
	analysis_freq_center_18 = 4350
	analysis_freq_center_19 = 5000
	analysis_freq_center_20 = 5675
	analysis_freq_center_21 = 6500
	analysis_freq_center_22 = 7500
	

	analysis_freq_high_1  = 313
	analysis_freq_high_2  = 438
	analysis_freq_high_3  = 563
	analysis_freq_high_4  = 688
	analysis_freq_high_5  = 813
	analysis_freq_high_6  = 938
	analysis_freq_high_7  = 1063
	analysis_freq_high_8  = 1188
	analysis_freq_high_9  = 1313
	analysis_freq_high_10 = 1563
	analysis_freq_high_11 = 1813
	analysis_freq_high_12 = 2063
	analysis_freq_high_13 = 2313
	analysis_freq_high_14 = 2688
	analysis_freq_high_15 = 3063
	analysis_freq_high_16 = 3563
	analysis_freq_high_17 = 4063
	analysis_freq_high_18 = 4688
	analysis_freq_high_19 = 5313
	analysis_freq_high_20 = 6063
	analysis_freq_high_21 = 6938
	analysis_freq_high_22 = 7938

	for thisChannel from 1 to numberOfChannels
	   # get cochlear position of center of analysis channel
		frequency = analysis_freq_center_'thisChannel'
		analysis_pos_center_'thisChannel' = log10((frequency/'aA')+'k')*'length'/'a'
		
		centerPos'thisChannel' = log10((analysis_freq_center_'thisChannel'/'aA')+'k')*'length'/'a'

	   # set electrode position (carrier) in the frequency domain, 
	   # including place shift
		carrier_freq_center_'thisChannel' = 'aA'*((10^('a'*(analysis_pos_center_'thisChannel'+shift)/'length'))-'k')
	endfor
endproc

procedure simulate_AB_ptp
   # match channel corner frequencies for the AB device
   # when using current focusing (partial tripolar mode)
   # as indicated on the experimental mapping software. 
	numberOfChannels = 14
	numberStimulated = 14
	
	analysis_freq_low_1  = 238
	analysis_freq_low_2  = 442
	analysis_freq_low_3  = 578
	analysis_freq_low_4  = 646
	analysis_freq_low_5  = 782
	analysis_freq_low_6  = 986
	analysis_freq_low_7  = 1189
	analysis_freq_low_8  = 1393
	analysis_freq_low_9  = 1665
	analysis_freq_low_10 = 2005
	analysis_freq_low_11 = 2413
	analysis_freq_low_12 = 2889
	analysis_freq_low_13 = 3500
	analysis_freq_low_14 = 4180

	analysis_freq_center_1  = 329
	analysis_freq_center_2  = 506
	analysis_freq_center_3  = 611
	analysis_freq_center_4  = 711
	analysis_freq_center_5  = 879
	analysis_freq_center_6  = 1083
	analysis_freq_center_7  = 1287
	analysis_freq_center_8  = 1523
	analysis_freq_center_9  = 1828
	analysis_freq_center_10 = 2200
	analysis_freq_center_11 = 2641
	analysis_freq_center_12 = 3180
	analysis_freq_center_13 = 3825
	analysis_freq_center_14 = 5810

	analysis_freq_high_1  = 442
	analysis_freq_high_2  = 578
	analysis_freq_high_3  = 646
	analysis_freq_high_4  = 782
	analysis_freq_high_5  = 986
	analysis_freq_high_6  = 1189
	analysis_freq_high_7  = 1393
	analysis_freq_high_8  = 1665
	analysis_freq_high_9  = 2005
	analysis_freq_high_10 = 2413
	analysis_freq_high_11 = 2889
	analysis_freq_high_12 = 3500
	analysis_freq_high_13 = 4180
	analysis_freq_high_14 = 8054

	for thisChannel from 1 to numberOfChannels
	   # get cochlear position of center of analysis channel
		frequency = analysis_freq_center_'thisChannel'
		analysis_pos_center_'thisChannel' = log10((frequency/'aA')+'k')*'length'/'a'
		
		centerPos'thisChannel' = log10((analysis_freq_center_'thisChannel'/'aA')+'k')*'length'/'a'

	   # set electrode position (carrier) in the frequency domain, 
	   # including place shift
		carrier_freq_center_'thisChannel' = 'aA'*((10^('a'*(analysis_pos_center_'thisChannel'+shift)/'length'))-'k')
	endfor
endproc


procedure simulate_custom_frequencies
   # Here's where you can fill in your own custom map
   # you can be as wild as you want. 
	numberOfChannels = 15
	numberStimulated = 15
	
	analysis_freq_low_1  = 250
	analysis_freq_low_2  = 421
	analysis_freq_low_3  = 505
	analysis_freq_low_4  = 607
	analysis_freq_low_5  = 730
	analysis_freq_low_6  = 877
	analysis_freq_low_7  = 1053
	analysis_freq_low_8  = 1266
	analysis_freq_low_9  = 1521
	analysis_freq_low_10 = 1827
	analysis_freq_low_11 = 2196
	analysis_freq_low_12 = 2638
	analysis_freq_low_13 = 3170
	analysis_freq_low_14 = 3809
	analysis_freq_low_15 = 4577

	analysis_freq_center_1  = 336
	analysis_freq_center_2  = 463
	analysis_freq_center_3  = 556
	analysis_freq_center_4  = 668
	analysis_freq_center_5  = 804
	analysis_freq_center_6  = 965
	analysis_freq_center_7  = 1160
	analysis_freq_center_8  = 1394
	analysis_freq_center_9  = 1674
	analysis_freq_center_10 = 2012
	analysis_freq_center_11 = 2417
	analysis_freq_center_12 = 2904
	analysis_freq_center_13 = 3490
	analysis_freq_center_14 = 4193
	analysis_freq_center_15 = 6638

	analysis_freq_high_1  = 421
	analysis_freq_high_2  = 505
	analysis_freq_high_3  = 607
	analysis_freq_high_4  = 730
	analysis_freq_high_5  = 877
	analysis_freq_high_6  = 1053
	analysis_freq_high_7  = 1266
	analysis_freq_high_8  = 1521
	analysis_freq_high_9  = 1827
	analysis_freq_high_10 = 2196
	analysis_freq_high_11 = 2638
	analysis_freq_high_12 = 3170
	analysis_freq_high_13 = 3809
	analysis_freq_high_14 = 4577
	analysis_freq_high_15 = 8700

	for thisChannel from 1 to numberOfChannels
	   # get cochlear position of center of analysis channel
		frequency = analysis_freq_center_'thisChannel'
		analysis_pos_center_'thisChannel' = log10((frequency/'aA')+'k')*'length'/'a'
		
		centerPos'thisChannel' = log10((analysis_freq_center_'thisChannel'/'aA')+'k')*'length'/'a'

	   # set electrode position (carrier) in the frequency domain, 
	   # including place shift
		carrier_freq_center_'thisChannel' = 'aA'*((10^('a'*(analysis_pos_center_'thisChannel'+shift)/'length'))-'k')
	endfor
endproc

procedure setParameters
   # convert names from the opening popup form, etc. 
   
   # temporal envelope filter
      envelope = envelope_cutoff_filter
   
   # basal shift
      shift = basal_shift
   
   # quick channel number sanity check
      if numberStimulated > numberOfChannels
  	 exit Number of stimulated channels cannot be greater than the number of total analysis channels!
      endif

   # cochlea parameters based on Greenwood (1990)
	aA = 165.4
	a = 2.1
	length = 35
	k = 0.88

   # potentially use this to set a minimum current width,
   # in the case of simulating output of a cochlear implant 
   # it's currently not used in any part of this script, 
   # but maybe something to use for the future. 
   	electrodeWidth = 0.75


 # ensure that window length for n-of-m processing is long enough
 # to capture the minimum pitch set by the envelope filter
 #interval_required = 6.4 / envelope
 #if interval < interval_required
 #   interval = interval_required
 #   print Changed time window to 'interval_required:3' to accomodate envelope filter'newline$'
 # endif

remove_vocoded_channel = 1
remove_sound_channel = 1
       
   # pre-emphasis and de-emphasis to aid in peak-picking
   # this accounts for the fact that speech energy rolls off 
   # by roughly 6 dB/octave,
   # resulting in disproportionate selection of low-frequency channels
   # in peak-picking mode. 
   # by preemphasizing, you're increasing spectral tilt by +6 dB/octave,
   # meaning the whole speech range shoudl be eligible for peak-picking. 
  	control.for.tilt = 1
       
   # re-name the rolloff parameter so that its string representation is fixed width
      if rolloff.per.mm < 10
   	rolloff.per.mm$ = "0'rolloff.per.mm'"
        else 
     	rolloff.per.mm$ = "'rolloff.per.mm'"
      endif
      
# decide on specific custom parameters (see end of script)
# if none chosen, work from basic Greenwood map
if filter = 4
	call simulate_custom_frequencies
   elsif filter = 3
   	call simulate_cochlear
   elsif filter = 2
   	call simulate_AB_ptp
   else
      # If you didn't select a pre-made custom map,
      # create channel frequency allocation
      # using channel numbers and freqeuncy ranges
      # typed in the startup window here.
	call setChannelCornerFrequencies lowCornerFreq highCornerFreq numberOfChannels
endif

if shape == 1
   shape$ = "peaked"
elsif shape == 2
   shape$ = "square"
elsif shape == 3
   shape$ = "LNN"
endif

if carrier_type == 1
    carrier_type$ = "noise"

elsif carrier_type == 2
   if shape$ != "LNN"
	   carrier_type$ = "toneComplex"
	    if pulseRate > 300
	       exit pulse rate too high! please choose a rate at 300 Hz or below
	    endif
    else
	   carrier_type$ = "noise"
	   print Changed carrier type to "noise" for LNN carriers'newline$'
    endif

elsif carrier_type ==3
    carrier_type$ = "sinewave"
    # override shape
    shape$ = "sine"

elsif carrier_type ==4
    carrier_type$ = "PSHC"
    print Pulse rate 'pulseRate' pps'newline$'

elsif carrier_type ==4
    carrier_type$ = "original_sound"
    print using original sound to create synthesis channels'newline$'
endif

if scramble_mix < 0 
	exit scramble_mix must be between 0 and 1 !
endif
if scramble_mix > 1
	exit scramble_mix must be between 0 and 1 !
endif

#**********************************************#
# create a suffix to the end of the object name 
# so you know which the vocoder settings
   if shape == 1
   	shape_label$ = "'rolloff.per.mm$'_dBmm"
   elsif shape == 2
   	shape_label$ = "flat_chan"
   elsif shape == 3
   	shape_label$ = "LNN"
   endif
   if carrier_type$ == "sinewave"
   # override that if it's a sinewave
   	shape_label$ = "sine"
   endif
   
   vocodedSuffix$ = "_voc_'numberStimulated'_of_'numberOfChannels'ch_'shape_label$'"
   
   # remove individual channels
   # (leave at 1 unless you're doing a demonstration with a single object
   # to create individual extracted channels)
   removeChannels = 1
   removeVocodedChannels = 1
      
   # declare pulse rate in final object name
   if carrier_type$ = "toneComplex"
   	vocodedSuffix$ = vocodedSuffix$ + "_pulse_'pulseRate:0'"
   endif
   
   keep_these_channels$ = "all"
   
   if keep_these_channels$ = "odd"
      	vocodedSuffix$ = vocodedSuffix$ + "_'keep_these_channels$'"
    elsif keep_these_channels$ = "even"
      	vocodedSuffix$ = vocodedSuffix$ + "_'keep_these_channels$'"
   endif
   
   
      # add DETAILED parameter info to the vocoded suffix
      # add ENV cutoff to vocodedSuffix$
      	vocodedSuffix$ = vocodedSuffix$ + "_e_'envelope'"
      	   if envelope == 0
      	   	vocodedSuffix$ = vocodedSuffix$ + "_e_scramble"
      	   endif

   compression_mult_for_label = compression_mult*100
   compression_mult_for_label = 'compression_mult_for_label:0'
   compression_mult_label$ = "'compression_mult_for_label'"
   vocodedSuffix$ = vocodedSuffix$ + "_comp_'compression_mult_label$'"
   
   
      # add SHIFT to vocodedSuffix$
           if shift > 0
         	   shift_label_val = shift * 10
         	   shift_label_val$ = "'shift_label_val'"
      	        vocodedSuffix$ = vocodedSuffix$ + "_s_'shift_label_val$'"
      	   else
      	        # if no shift, just leave it blank
      	        vocodedSuffix$ = vocodedSuffix$
      	        # verbose suffix instead 
      	        vocodedSuffix$ = vocodedSuffix$ + "_s_0"
      	endif
      	
   # if you're scrambling the envleope, add that to the suffix as well
   #if envelope == 0 
   #	vocodedSuffix$ = vocodedSuffix$ + "_scrambleEnv"
   #endif
   
   # NOTE:
   # if you do not like those naming schemes, you can create a new one here, 
   # and just un-comment the line
   # vocodedSuffix$ = "your_vocoder_suffix_here"
   
   if remove_suffix == 1
   	vocodedSuffix$ = ""
   endif
   
   # hilbert envelope might slow you down a lot. 
   use_Hilbert = 0
   
   num_LNN_iterations = 4
   
   # print out detailed progress?
   verbose = 1
endproc