center_freq = 1000
duration = 3
num_iterations = 5
remove_iterations = 1
samplerate = 48000
stim_name$ = "my_LNN_1000"



clearinfo


call make_LNN_ERB center_freq duration num_iterations samplerate remove_iterations 'stim_name$'

procedure make_LNN_ERB cf .duration .num_iterations .samplerate .remove_iterations .newname$
  use_Hilbert = 1
  # print status of each iteration
  verbose = 1
  
  # 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_chunks 'target_sound$'
  		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
	
	if .remove_iterations == 1
	   for n from 1 to .num_iterations
	      index_to_remove = n - 1
	      select Sound 'original_sound$'_LNN_'index_to_remove'
	      Remove
	   endfor
	else
	   # sequnce the iterations
		selectObject: "Sound noise_to_filter_LNN_0"
		   for iteration_index from 1 to (.num_iterations)
			plusObject: "Sound noise_to_filter_LNN_'iteration_index'"
		   endfor
		Concatenate recoverably
		selectObject: "Sound chain"
		Rename: "LNN_iterations"
	   	selectObject: "TextGrid chain"
		Rename: "LNN_iterations"	
	
	endif
	

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

endproc



procedure get_envelope_TFS_chunks .name$

	# Get duration to be used later
	select Sound '.name$'
	.duration = Get total duration
	.samplerate = Get sampling frequency

	# divide into 32768-sampel chunks 
	# to facilitate faster Hilbert transform
	call divide_32768 '.name$'

	# figure out how many chunks were created
	num_chunks = chunks_analyzed
	
	# Hilbert transform for each successive chunk
	for sound_part from 1 to num_chunks
		call get_envelope_TFS '.name$'_part_'sound_part'
	endfor
	
	# Sequence the TFSs back together into one sound
	   nocheck select junk
	   for sound_part from 1 to num_chunks
	   	plus Sound '.name$'_part_'sound_part'_TFS
	   endfor
	   Override sampling frequency: .samplerate
	   Concatenate
	   Rename... junk
	   Extract part: 0, .duration, "rectangular", 1, "no"
	   Rename... '.name$'_TFS
	   select Sound junk
	   Remove
	
	# Sequence the ENVs back together into one sound
		nocheck select junk
		for sound_part from 1 to num_chunks
		   	plus Sound '.name$'_part_'sound_part'_ENV
		endfor
		Override sampling frequency: .samplerate
		Concatenate
		Rename... junk
		Extract part: 0, .duration, "rectangular", 1, "no"
		Rename... '.name$'_ENV
		select Sound junk
	   	Remove
	
	# cleanup
	nocheck select junk
	  for sound_part from 1 to num_chunks
	   	plus Sound '.name$'_part_'sound_part'
	   	plus Sound '.name$'_part_'sound_part'_TFS
		plus Sound '.name$'_part_'sound_part'_ENV
	  endfor
	  Remove

endproc



procedure get_envelope_TFS .name$
.lpf = 8000
  
    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 divide_32768 .name$

	select Sound '.name$'
	samplerate = Get sampling frequency
	duration = Get total duration
	num_samples = Get number of samples

	samples_collected = 0
	chunks_analyzed = 0
	sample_start = 0
	while samples_collected < num_samples

	sample_end = sample_start + 32768
	time_start = Get time from sample number: sample_start
	time_end = Get time from sample number: sample_end

	selectObject: "Sound '.name$'"
	Extract part: time_start, time_end, "rectangular", 1, "no"

	samples_collected = samples_collected + 32768
	chunks_analyzed = chunks_analyzed + 1
	sample_start = sample_start + 32769

	#print chunk 'chunks_analyzed' 'time_start:3' 'time_end:3''newline$'

	Rename... '.name$'_part_'chunks_analyzed'

	endwhile

endproc