# F0-controlled bandfilter analysis script for the program 'Praat'
# One F0 period of a voiced speech sound must be isolated and selected in advance;
# 	the exact period crossing need not be determined
# Unvoiced intervals can be chosen freely
# The output is an approximation of the envelope spectrum

form F0-controlled Gauss BF envelope spectrum
   comment Method: moving Gauss window on DFT components of sustained sound.
   positive Effective_filter_width_/F0_(Hz) 1.1
   positive Highest_frequency_(Hz) 6000
   boolean Pre-emphasis_6dB/oct no
   boolean Draw_spectrum yes
endform


# shorten variables
bwcorr = 'Effective_filter_width_/F0'
fhi = highest_frequency
pe = 'Pre-emphasis_6dB/oct'

obj$ = selected$ ("Sound")
d = Get duration

# check frequency range
f0 = 1/d
neededRange = fhi + f0/2
sfreq = Get sampling frequency
maxPossible_f = sfreq/2 - f0/2
if fhi > maxPossible_f
   exit The specified maximum frequency is too high.
   ... The sample frequency allows 'maxPossible_f' Hz only.
endif

# make sustained sound by concatenating periods
ns = Get number of samples
length = 3*d
Create Sound... sustained 0 'length' 'sfreq' Sound_'obj$'(x mod 'd')
Rename... segment

if pe = 1
     Pre-emphasize (in-line)... 50
endif

Multiply by window... Hanning
To Spectrum... no
Copy... mult

# calc. parameters for Gauss window
fbw = bwcorr * f0
asquare = 2*ln(2)/fbw^2

# set up Matrix for spectral values
fstep = fbw/3
numsteps = ceiling(fhi/fstep)+1
Create Matrix... 'obj$' 0 ('numsteps'-1)*'fstep' 'numsteps'
... 'fstep' 0 1 2 2 1 1 0

for step to numsteps
   f = (step-1) * fstep
   select Spectrum mult
   Formula... Spectrum_segment()*exp(-'asquare'*(x-'f')^2)
   energy = Get band energy... 0 0
   specComp = sqrt(energy*fbw)
   # filling one row suffices: second remains 0
   select Matrix 'obj$'
   Set value... 1 'step' 'specComp'
endfor

To Spectrum
if draw_spectrum = 1
   Draw... 0 'fhi' 0 0 yes
endif

# cleanup
select Sound segment
plus Spectrum segment
plus Spectrum mult
plus Matrix 'obj$'
Remove

# end of program