form Truncated-filter spectrum
     real Filter_width_/_F0 0.25
     positive Maximum_frequency_(Hz) 5000
endform


obj$ = selected$("Sound")
Copy... period

# Reverse to get cross correlation using convolution
Reverse

d = Get total duration
bw = 'Filter_width_/_F0'/d
fsamp = Get sample rate

# apply 5 steps per F0
fstep = 1/d/5

nfreq = ceiling(maximum_frequency/fstep)
Create Matrix... spectrum 0 maximum_frequency
... nfreq fstep 0 1 2 2 1 1 0

# calculate phase step
step1500 = 1500/fstep
phasestep = pi/2/step1500

for ifreq to nfreq
     f = (ifreq-1)*fstep
     if ifreq < step1500
          phase = (ifreq-1)*phasestep
     else
          phase = pi/2
     endif

     # create filter i.r.
     Create Sound from formula... ir Mono 0 'd' 'fsamp' sin(2*pi*'f'*x+'phase')*exp(-pi*'bw'*x)
     plus Sound period
     Convolve

     energy = Get energy... 0 'd'

     # remove temp objects
     select Sound ir
     plus Sound period_ir
     Remove

     #Store the result
     select Matrix spectrum
     Set value... 1 ifreq sqrt(energy)

endfor

To Spectrum

# calibrate
ampF0 = Get real value in bin... 5+1
select Sound period
To Spectrum... no
re = Get real value in bin... 2
im = Get imaginary value in bin... 2
amp = sqrt(re^2+im^2)
corr = amp/ampF0
Remove
select Spectrum spectrum
Formula... self * 'corr'

# smooth with moving average filter and correct for frequency shift
Formula... (self[col-1]+self[col]+self[col+1])/3

Rename... 'obj$'