Trying to get a csound trumpet synth from the 1995 journal of acoustics working

Newbies
#1

I’ve got a journal paper from the 1995 journal of acoustics with a listing for a wavetable trumpet synth. It’s in scanned pdf format, so have to type it in. What I’ve typed is the listing below.
Also, the definition of the kbr uses a function called “as”. I can’t find any reference to what that means. I wonder if that is supposed to be the “abs” function.

The program runs, but it crackles. I may have made some typing errors and there may be typing errors in the actual paper itself. I think I have already found a couple of errors.

Any pointers would be greatly appreciated.

-o dac // real-time output sr = 22050 // sample rate kr=2205 0dbfs = 1 // maximum amplitude (0 dB) is 1 nchnls = 2 // number of channels is 2 (stereo) ksmps = 10 // number of samples in one control cycle (audio vector) giseed=5.5

instr 2
// sine tone with amplitude=0.2 and frequency=550Hz
aSine = poscil:a(0.2,880)
// apply fade-in 0.1 sec and fade-out of 1 sec over the duration (p3)
aOut = linen:a(aSine,0.1,p3,1)
outall(aOut)
endin

instr 1
ipchfreq=p4
iamp1=p5
iamp2=p6
imaxbr=p7

ifreq=cpspch(ipchfreq)
iattack = 0.24*(1.1667 - .1667*ifreq/261)
if iattack<0.30*p3 igoto init1
    iattach=0.30*p3

init1: 
    idecay=0.8
    if idecay<0.65*p3 igoto init2
    idecay=0.65*p3  
       
init2: isustain=p3-iattack-idecay
    if isustain !=0 igoto init3
    isustain=.1*idecay
    idecay=.9*idecay
    
init3:
        kpct line 0.0, p3, 1.0
        ktime = kpct*p3
        kpt = kpct*100

;prototype dfr curve (from ctp03.pv.an)
    kfreq linseg-.015, iattack/2,-.004,
    iattack/2, .0, isustain, .0, idecay/2, .005,    
idecay/2, .021
   kfreqr1 randi .0004, 10, giseed         

giseed=frac(giseed+3.141592654)
kfreqr2 rand .002, giseed


giseed=frac(giseed+3.141592654)
kfreqr2 rand .002, giseed

giseed=frac(giseed+3.141592654)
kfreq=(kfreq+kfreqr1+kfreqr2)*ifreq
display kfreq, p3

;prototype rms curve (from ctpt03.pv.an)
krms linseg 0,.031*iattack,.021*iamp1,
.031*iattack, .028*iamp1, .041*iattack, .21*iamp1,
.01*iattack, .18*iamp1, .017*iattack, .28*iamp1,
.08*iattack, .4*iamp1, .16*iattack, .44*iamp1,
.32*iattack, .7*iamp1, .06*iattack, .714*iamp1,
.25*iattack, iamp1, isustain, iamp2, .26*idecay,
.075*iamp2, .18*idecay, .42*iamp2, .31*idecay,
.109*iamp2, .17*idecay, .0109*iamp2, .015*idecay,
.0028*iamp2, .03*idecay, .0019*iamp2, .005*idecay,
0
    krmsr randi .05*krms, 20, giseed
    
    giseed =frac(giseed+3.141592654)
    krms = abs(krms+krmsr)
    display krms, p3
    
    kbr linseg 0, .01*iattack, .33*imaxbr,
    .1*iattack, .48*imaxbr, .3*iattack,
    .68*imaxbr, .4*iattack, .96*imaxbr, .1*iattack+
    .1*isustain, 1.0*imaxbr, .6*isustain, .94*imaxbr,
    .3*isustain, .54*imaxbr, .21*idecay, .28*imaxbr,
    .16*idecay, .15*imaxbr, .16*idecay, 0
    
    kbrr1 randi 0.07*kbr, 20, giseed
    
    giseed=frac(giseed+3.141592654)
    kbrr2 rand 0.02*kbr, giseed
    
    kbrr2 = 0
       giseed=frac(giseed+3.141592654)
       //kbr=as(kbr+kbrr1+kbrr2)
       kbr = kbrr1+kbrr2
       display kbr, p3
       
    w1: if ipchfreq>8.01 goto w2
        iwt = 1
        goto wend
    w2: if ipchfreq>8.0 goto w3
        iwt = 2
        goto wend
    w3: if ipchfreq>9.01 goto w4
        iwt = 3
        goto wend
    w4: if ipchfreq>10.00 goto w5
        iwt = 4
        goto wend
    w5: iwt = 5
    wend: awt oscili krms, ifreq+kfreq, iwt, giseed
          giseed=frac(giseed+3.141592654)
          
    ; set filter parameters
        icut0 = 20
        icut1 = 300
        icut2 = 500
        icut3 = 825
        icut4 = 975
        icut5 = 1225
        icut6 = 1525
        icut7 = 1925  
        icut8 = 2450
        icut9 = 3000
        icut10 = 4000
        
        kw2 = frac(kbr)
        kw1 = 1.0-kw2
     f1: if kbr>=1 kgoto f2
            kcut = kw1*icut0+ kw2*icut1
            kgoto fend
     f2: if kbr>=2 kgoto f3
            kcut = kw1*icut1 + kw2*icut2
            kgoto fend
     f3: if kbr>=3 kgoto f4
            kcut = kw1*icut2 + kw2*icut3
            kgoto fend
     f4: if kbr>=4 kgoto f5
            kcut = kw1*icut3 + kw2*icut4
            kgoto fend
     f5: if kbr>=5 kgoto f6
            kcut = kw1*icut4 + kw2*icut5
            kgoto fend
     f6: if kbr>=6 kgoto f7
            kcut = kw1*icut5 + kw2*icut6
            kgoto fend
     f7: if kbr>=7 kgoto f8
            kcut = kw1*icut6 + kw2*icut7
            kgoto fend
     f8: if kbr>=8 kgoto f9
            kcut = kw1*icut7 + kw2*icut8
            kgoto fend
     f9: if kbr>=9 kgoto f10
            kcut = kw1*icut8 + kw2*icut9
            kgoto fend
     f10: if kbr>=10 kgoto fmax
             kcut = kw1*icut9 + kw2*icut10
             kgoto fend
     fmax: kcut = icut10
     fend:
             display kcut, p3
             afilt reson awt, 0, sqrt(2.)*kcut, 0
             asig balance afilt, awt
             kgoto end
     end:
             anoise rand 0.001*asig, giseed
             giseed = frac(giseed+3.141592654)
             asig = asig+anoise
             out asig

endin

// call instrument 1 at start time 0 and for 2 seconds i 2 0 2

f1 0 8193 10 0.376492 -0.836273 0.875826
-0.702022 -0.944316 0.782612 0.755986 0.755986
0.785038 0.588679 -0.544976 0.544976 0.489619
0.489619 -0.233817 -233817 0.132038
-0.132038 0.132038 -0.072570 0.072570 0.072570
0.114808 0.114808 -0.114808 -0.114808 0.059510
0.059510 0.059510 0.059510 -0.059510 0.059510
0.049870 0.049870 0.049870 -0.049870 0.049870
-0.049870 0.049870 0.01302 0.013012 0.013012

f2 0 8193 10 0.574816 -0.875826 -0.702022
0.782712 -0.755986 -0.785038 -0.588679
-0.544976 0.489619 -0.233817 0.233817
-0.132038 -0.072570 0.072570 -0.072570
-0.114808 -0.114808 -0.114808 0.059510
-0.059510 -0.059510 0.049870 0.049870 0.049870
-0.049870 -0.049870 0.013012 0.013012 0.013012

f3 0 8193 10 08836273 0.702022 0.782612 -0.785038
0.588679 0.544976 -0.489619 0.233817 0.132038
0.072570 -0.072570 -0.114808 -0.114808
0.059510 -0.59510 0.049870 0.049870 -0.049870
-0.049870 -0.013012 -0.013012

f4 0 8193 10 0.875826 0.782612 -0.785038 0.544976
-0.233817 -0.132038 -0.072570 -0.114808
-0.114808 0.059510

f5 0 8193 10 0.875826 0.782612 -0.544976 0.233817
-0.132038 -0.072570 -0.022570

i1 0 3.0 8.01 2000 1800 9
i1 2 3.0 7.08 2000 1800 9
i1 2 3.0 9.01 2000 1800 9
i1 4 3.0 7.05 2000 1800 9
i1 4 3.0 8.08 2000 1800 9
i1 4 3.0 9.00 2000 1800 9

#2

Hi George,

I found the article on the web and checked the original code. It contains some typos but I think I corrected them well. I also added my comments on the unclear points which I can’t solve.
The great thing about Csound is that the code from 30 years ago works as is without any changes!

<CsoundSynthesizer>
<CsLicense>
Synthesis of Trumpet Tones Using a Wavetable and a Dynamic filter
Andrew Horner and James Beauchamp
J. Audio Eng. Soc., Vol. 43, No. 10, 1995 October
https://cmp.ischool.illinois.edu/publications.html
</CsLicense>
<CsOptions>
-odac
</CsOptions>
<CsInstruments>
sr = 22050 // sample rate
kr = 2205
ksmps = 10 // number of samples in one control cycle (audio vector)
nchnls = 2 // number of channels is 2 (stereo)

;0dbfs = 32767
;The article is from 1995, which is before the release of Csound v4.10 when odbfs was introduced. Accordingly 0dbfs should be 32767, but can be omitted.

giseed=5.5

instr 1
ipchfreq=p4
iamp1=p5
iamp2=p6
imaxbr=p7

ifreq=cpspch(ipchfreq)
iattack = 0.24*(1.1667 - .1667*ifreq/261)
;don't allow attack to be more than 30% of tone
if iattack<0.30*p3 igoto init1
    iattack=0.30*p3

init1:
    idecay=0.8
;don't allow decay to be more than 65% of tone
    if idecay<0.65*p3 igoto init2
    idecay=0.65*p3 

init2:
    isustain=p3-iattack-idecay    ;typo?(decay -> idacay)
    if isustain !=0 igoto init3
    isustain=.1*idecay
    idecay=.9*idecay

init3:
    kpct line 0.0, p3, 1.0
    ktime = kpct*p3
    kpt = kpct*100
    ;Note: kpct,ktime,kpt are not used anywhere

;prototype dfr curve (from ctp03.pv.an)
    kfreq linseg -.015, iattack/2,-.004, iattack/2, .0, isustain, .0, idecay/2, .005, idecay/2, .021
    kfreqr1 randi .0004, 10, giseed
;0.4% nuances
    giseed=frac(giseed+3.141592654)
    kfreqr2 rand .002, giseed
;0.2% nuanes
    giseed=frac(giseed+3.141592654)
    kfreq=(kfreq+kfreqr1+kfreqr2)*ifreq
    display kfreq, p3
;prototype rms curve (from ctpt03.pv.an)
    krms linseg 0,.031*iattack,.021*iamp1,.031*iattack, .028*iamp1, .041*iattack, .21*iamp1,.01*iattack, .18*iamp1, .017*iattack, .28*iamp1,.08*iattack, .4*iamp1, .16*iattack, .44*iamp1,.32*iattack, .7*iamp1, .06*iattack, .714*iamp1,.25*iattack, iamp1, isustain, iamp2, .26*idecay,.75*iamp2, .18*idecay, .42*iamp2, .31*idecay,.109*iamp2, .17*idecay, .0109*iamp2, .015*idecay,.007*iamp2,.01*idecay,.0047*iamp2,.02*idecay,.0028*iamp2,.03*idecay,.0019*iamp2,.005*idecay,0
    krmsr randi .05*krms, 20, giseed
    ;Note: I guess that the 25th parameter should be ".75*iamp2" though I don't check the details of the article.
    ;The original code ".075*iamp2)" results in non-contiguous envelope.
    
;5% nuances
    giseed =frac(giseed+3.141592654)
    ;krms = abs(rksm+krmsr)   ;typo?(rksm ->krms)
    krms = abs(krms+krmsr)
    display krms, p3

;prototype pseudo-brightness curve
    kbr linseg 0, .01*iattack, .33*imaxbr, .1*iattack, .48*imaxbr, .3*iattack, .68*imaxbr, .4*iattack, .96*imaxbr, .1*iattack+ .1*isustain, 1.0*imaxbr, .6*isustain, .94*imaxbr, .3*isustain, .85*imaxbr,.26*idecay, .70*imaxbr,.21*idecay,.54*imaxbr, .21*idecay, .28*imaxbr, .16*idecay, .15*imaxbr, .16*idecay, 0
    kbrr1 randi 0.07*kbr, 20, giseed
;7% nuances
    giseed=frac(giseed+3.141592654)
    kbrr2 rand 0.02*kbr, giseed
    ;Note: kbrr2 is always 0 because it is overwritten in the next row
;2% noise
    kbrr2 = 0
    giseed=frac(giseed+3.141592654)
    ;kbr=as(kbr+kbrr1+kbrr2)   ;typo?(as -> abs)
    kbr=abs(kbr+kbrr1+kbrr2)
    display kbr, p3
    w1: if ipchfreq>8.01 goto w2
        iwt = 1
        goto wend
    w2: if ipchfreq>8.0 goto w3
        iwt = 2
        goto wend
    w3: if ipchfreq>9.01 goto w4
        iwt = 3
        goto wend
    w4: if ipchfreq>10.00 goto w5
        iwt = 4
        goto wend
    w5: iwt = 5
    wend: awt oscili krms, ifreq+kfreq, iwt, giseed
          giseed=frac(giseed+3.141592654)

; set filter parameters
       icut0 = 20
       icut1 = 300
       icut2 = 500
       icut3 = 825
       icut4 = 975
       icut5 = 1225
       icut6 = 1525
       icut7 = 1925
       icut8 = 2450
       icut9 = 3000
       icut10 = 4000

       kw2 = frac(kbr)
       kw1 = 1.0-kw2
f1: if kbr>=1 kgoto f2
       kcut = kw1*icut0+ kw2*icut1
       kgoto fend
f2: if kbr>=2 kgoto f3
       kcut = kw1*icut1 + kw2*icut2
       kgoto fend
f3: if kbr>=3 kgoto f4
       kcut = kw1*icut2 + kw2*icut3
       kgoto fend
f4: if kbr>=4 kgoto f5
       kcut = kw1*icut3 + kw2*icut4
       kgoto fend
f5: if kbr>=5 kgoto f6
       kcut = kw1*icut4 + kw2*icut5
       kgoto fend
f6: if kbr>=6 kgoto f7
       kcut = kw1*icut5 + kw2*icut6
       kgoto fend
f7: if kbr>=7 kgoto f8
       kcut = kw1*icut6 + kw2*icut7
       kgoto fend
f8: if kbr>=8 kgoto f9
       kcut = kw1*icut7 + kw2*icut8
       kgoto fend
f9: if kbr>=9 kgoto f10
       kcut = kw1*icut8 + kw2*icut9
       kgoto fend
f10: if kbr>=10 kgoto fmax
       kcut = kw1*icut9 + kw2*icut10
       kgoto fend
fmax: kcut = icut10
fend:
       display kcut, p3
       afilt reson awt, 0, sqrt(2.)*kcut, 0
       asig balance afilt, awt
       kgoto end
end:
;.1% noise in signal
       anoise rand 0.001*asig, giseed
       giseed = frac(giseed+3.141592654)
       asig = asig+anoise
       out asig
endin
</CsInstruments>
<CsScore>
// call instrument 1 at start time 0 and for 2 seconds i 2 0 2
f1 0 8193 10 0.376492 -0.836273 0.875826 -0.702022 -0.944316 0.782612 0.755986 0.755986 0.785038 0.588679 -0.544976 0.544976 0.489619 0.489619 -0.233817 -0.233817 0.132038 -0.132038 0.132038 -0.072570 0.072570 0.072570 0.114808 0.114808 -0.114808 -0.114808 0.059510 0.059510 0.059510 0.059510 -0.059510 0.059510 0.049870 0.049870 0.049870 -0.049870 0.049870 -0.049870 0.049870 0.01302 0.013012 0.013012

f2 0 8193 10 0.574816 -0.875826 -0.702022 0.782612 -0.755986 -0.785038 -0.588679 -0.544976 0.489619 -0.233817 0.233817 -0.132038 -0.072570 0.072570 -0.072570 -0.114808 -0.114808 -0.114808 0.059510 -0.059510 -0.059510 0.049870 0.049870 0.049870 -0.049870 -0.049870 0.013012 0.013012 0.013012

f3 0 8193 10 0.836273 0.702022 0.782612 -0.785038 0.588679 0.544976 -0.489619 0.233817 0.132038 0.072570 -0.072570 -0.114808 -0.114808 0.059510 -0.59510 0.049870 0.049870 -0.049870 -0.049870 -0.013012 -0.013012

f4 0 8193 10 0.875826 0.782612 -0.785038 0.544976 -0.233817 -0.132038 -0.072570 -0.114808 -0.114808 0.059510

f5 0 8193 10 0.875826 0.782612 -0.544976 0.233817 -0.132038 -0.072570 -0.022570

i1 0 3.0 8.01 2000 1800 9
i1 2 3.0 7.08 2000 1800 9
i1 2 3.0 9.01 2000 1800 9
i1 4 3.0 7.05 2000 1800 9
i1 4 3.0 8.08 2000 1800 9
i1 4 3.0 9.00 2000 1800 9

</CsScore>
</CsoundSynthesizer>
#3

Thanks Tetsuya, that works. I did manage to get my attempt to run eventually, but I hadn’t fixed some of the errors that you fixed.
I’ll have a closer look at what you’ve done over the next few days.
Just wondering if you have any ideas of what the best way to display the brightness and amplitude envelopes? I am using Cabbage and they have some widgets for displaying waveforms, but I wanted to look at the shape of the envelopes in this trumpet simulation.

#4

You can output any k-sig to an audio file and open it with audio editor like Audacity.
I actually did it to check the envelope of krms.

For example, try changing the code as follows:
out asig
→ out a(krms)
And change CsOptions:
-odac
→ -W -o envelope_krms.wav

Please use CsoundQT or run the csd in terminal because rendering to audio file is not possible by Cabbage.