# rotation.praat
# creates one low-pass filtered version
# and one low-pass filtered and spectrally rotated version
# from all sound files in a given directory (<- should contain no other files)

# rotation around 0.5*MaximumFrequency
# low-pass filtering at 0.95*MaximumFrequency
# rms scaled to match the input sound
# ltas adjusted to approximate the input sound

# ltas adjustment written by Holger Mitterer
# and made faster with the help of Paul Boersma
# frank eisner 07/05


# clear object and info windows
select all
nocs = numberOfSelected ("Sound")
if nocs >= 1
  Remove
endif
clearinfo

### Input form
form Spectral Rotation
positive Maximum_frequency 4000
word directory K:\sounds\pilot_ABX\temp\
positive LTAS_bandwidth 100
boolean Write_files 0
endform

wd$=directory$
maxfreq = 'Maximum_frequency'
ltasbw = 'LTAS_bandwidth'

# # # main loop
printline 'wd$'
Create Strings as file list... list 'wd$'*.wav
nfiles = Get number of strings
for cmfile to nfiles
	select Strings list
	sound$ = Get string... cmfile
	Read from file... 'wd$''sound$'
	fname$ = selected$ ("Sound", 1)


printline Sound 'cmfile' of 'nfiles'

# rotate then adjust ltas
call rotation 'fname$' 'maxfreq'
Rename... Rotation only
Copy... 'fname$'_r'maxfreq'_band
call adjust_ltas 'fname$'_l'maxfreq' 'fname$'_r'maxfreq'_band 


# clean up

select Sound 'fname$'_r'maxfreq'
plus Sound org
plus Spectrum 'fname$'_r'maxfreq'_band
plus Ltas 'fname$'_r'maxfreq'_band
plus Spectrum 'fname$'_l'maxfreq'
plus Ltas 'fname$'_l'maxfreq'
plus Spectrum org
plus Spectrum toFilter
plus Spectrum temp
plus Ltas diffltas
Remove
select Sound 'fname$'_r'maxfreq'_band__av
Rename... 'fname$'_r'maxfreq'
Copy... adjusted 
call adjust_intensity_contour 'fname$' 'maxfreq'

# write sound files
if 'write_files' == 1
select Sound 'fname$'_r'maxfreq'
Write to WAV file... 'wd$''fname$'_r'maxfreq'.wav
select Sound 'fname$'_l'maxfreq'
Write to WAV file... 'wd$''fname$'_l'maxfreq'.wav
endif

printline

# # # close main loop
endfor

# clean up
select Strings list
Remove

# write to info window
if 'write_files' == 1
printline
printline Output files written to 'wd$'.
endif
printline
printline Done.

### End


# ====================================================================
# rotation
# creates one low-pass filtered version and one low-pass filtered 
# rotated version; scales both to match RMS intensity of the input file
# ====================================================================

procedure rotation rname$ maxfreq

select Sound 'rname$'
orms = Get root-mean-square... 0 0
sf = Get sample rate
hsf = sf/2
ttime = Get total duration
stopfreq = 'maxfreq' - 150

if maxfreq > hsf
	exit Maximum frequency must be less than 'hsf' Hz
endif

# Low-pass filter original and scale RMS
printline ...low-pass filtering
select Sound 'rname$'
Filter (pass Hann band)... 150 'stopfreq' 100
lrms = Get root-mean-square... 0 0
Formula... self * (orms/lrms)
Rename... 'rname$'_l'maxfreq'

# Rotate and scale RMS
printline ...rotating
select Sound 'rname$'_l'maxfreq'
Copy... 'rname$'_r'maxfreq'
Formula... self * 1/2 * sin(2*pi*'maxfreq'*x)
Filter (pass Hann band)... 0 'stopfreq' 100
rrms = Get root-mean-square... 0 0
Formula... self * (orms/rrms)

endproc


# ====================================================================
# adjust_ltas
# filters the rotated sound so that the long-term average spectrum
# is approximately equal to that of the input sound.
#====================================================================

procedure adjust_ltas insnd$ outsnd$

# playback final sound if set to 1
  listenIn = 0

  printline ...adjusting LTAS

# create LTAS of rotated and unrotated sounds
   select Sound 'outsnd$'
   To Spectrum (dft)
   To Ltas... 'ltasbw'
   select Sound 'insnd$'
   To Spectrum (dft)
   To Ltas... 'ltasbw'
# compute the difference
   Copy... diffltas
   Formula... self[col] - Ltas_'outsnd$'[col]


### to-be filtered sound
   select Sound 'outsnd$'
   t = Get total duration
   sr = Get sampling frequency
   st = Get starting time
   et = Get finishing time

   Rename... org

# generate its spectrum
   To Spectrum (dft)
   select Spectrum org    
   Copy... toFilter

### filter
  select Spectrum toFilter
  Formula... self * 10 ^ (Ltas_diffltas (x) / 20)
  Copy... temp
  To Sound (fft)
   Rename... core
# resample
   select Sound core
   Rename... temp
   Resample... 'sr' 50
   Rename... final
   select Sound temp
   Remove

   select Sound final
   Rename... 'outsnd$'_av
   if listenIn
	Play
   endif

endproc

# ====================================================================
# adjust_intensity_contour
#====================================================================

procedure adjust_intensity_contour speech$ 

  select Sound 'speech$'_l'maxfreq'

  # Get the intensity value in dB
  origIntensityValue = Get intensity (dB)
  To Intensity... 70 0
  Rename... origIntensityObj
  Down to IntensityTier
  Rename... origIntensityTierObj

  # Intensity object of the rotated utterance
  select Sound 'speech$'_r'maxfreq'
  To Intensity... 70 0
  Rename... rotatedIntensityObj1

  # Inverse the intensity object by getting the maximum and subtracting self
  maxSyn = Get maximum... 0 0 Parabolic
  Formula... 'maxSyn' - self

  # And make IntensityTier object
  Down to IntensityTier
  Rename... rotatedIntensityTierObj1

  # Multiply the synthetic utterance with its own inverse IntensityTier and then
  # by the IntensityTier of the natural utterance
  select Sound 'speech$'_r'maxfreq'
  plus IntensityTier rotatedIntensityTierObj1
  Multiply
  Rename... even
  plus IntensityTier origIntensityTierObj
  Multiply
  Rename... completeSoundObj1



#clean up
   select Intensity origIntensityObj
   plus IntensityTier origIntensityTierObj
   plus Intensity rotatedIntensityObj1
   plus IntensityTier rotatedIntensityTierObj1
   plus Sound even
   plus Sound 'speech$'_r'maxfreq'
   Remove
   select Sound completeSoundObj1
   Rename... 'speech$'_r'maxfreq'
 # Before writing, adjust the average intensity value in dB
 # add fade-in and -out and scale to match orig rms
    printline ...scaling    
    endtime = Get total duration
    Formula... if x < 0.005 then self * (x/0.005) else self fi
    Formula... if x > ('endtime') - 0.005 then self * (('endtime'-x)/0.005) else self fi      
    Scale intensity... 'origIntensityValue'
endproc