parcoe.for

C       PARCOE.FOR             D.H. KLATT               8/1/78
C
C       "PARAM-TO-COEF" TRANSFORMATION SUBROUTINE
C
C     THIS PROGRAM CONVERTS SYNTHESIZER CONTROL PARAMETERS FROM ARRAY I(39)
C     INTO DIFFERENCE EQUATION CONSTANTS FOR SYNTHESIZER HARDWARE
C     STORED IN ARRAY C(50)
C
        SUBROUTINE PARCOE(INITPC)
C               INITPC INITIALIZES THIS ROUTINE IF =-1

        use, intrinsic:: iso_fortran_env, only: int16, sp=>real32

        IMPLICIT INTEGER(int16) (I-N)
        IMPLICIT REAL(sp) (A-D, S, T, P)

        implicit none (external)

        integer(int16), intent(inout) :: initpc

        real(sp), external :: GETAMP
        external :: SETABC

        REAL(sp) :: IMPULS, C(50)
        INTEGER(int16) I(39),NDBSCA(12),NDBCOR(10)

        real(sp), parameter :: pi=4.*atan(1.)
C     INPUT PARAMETER VALUES (CONSTANT AND VARIABLE) PASSED THROUGH I
        COMMON /PARS/ I
        COMMON /COEFS/ C
        COMMON /PIXX/ PIT,TWOPIT
C     COEFICIENT VALUES IN C(50) ARE REAL
        EQUIVALENCE (C(1),IMPULS),(C(2),SINAMP),(C(3),AFF)
        EQUIVALENCE (C(4),AHH),(C(5),A1P),(C(6),A2P)
        EQUIVALENCE (C(7),A3P),(C(8),A4P),(C(9),A5P)
        EQUIVALENCE (C(10),A6P),(C(11),ABP),(C(12),ANPP)
        EQUIVALENCE (C(13),AGP),(C(14),BGP),(C(15),CGP)
        EQUIVALENCE (C(16),AGZ),(C(17),BGZ),(C(18),CGZ)
        EQUIVALENCE (C(19),AGS),(C(20),BGS),(C(21),CGS)
        EQUIVALENCE (C(22),A1),(C(23),B1),(C(24),C1)
        EQUIVALENCE (C(25),A2),(C(26),B2),(C(27),C2)
        EQUIVALENCE (C(28),A3),(C(29),B3),(C(30),C3)
        EQUIVALENCE (C(31),A4),(C(32),B4),(C(33),C4)
        EQUIVALENCE (C(34),A5),(C(35),B5),(C(36),C5)
        EQUIVALENCE (C(37),A6),(C(38),B6),(C(39),C6)
        EQUIVALENCE (C(40),ANP),(C(41),BNP),(C(42),CNP)
        EQUIVALENCE (C(43),ANZ),(C(44),BNZ),(C(45),CNZ)
        EQUIVALENCE (C(46),PLSTEP)
C     NAME OF INPUT CONTROL PARAMETERS
        EQUIVALENCE (I(1),NNAV),(I(2),NNAF),(I(3),NNAH),(I(4),NNAVS)
     1,(I(5),NNF0),(I(6),NNF1),(I(7),NNF2),(I(8),NNF3),(I(9),NNF4)
     1,(I(10),NNFNZ),(I(11),NNAN),(I(12),NNA1),(I(13),NNA2)
     1,(I(14),NNA3),(I(15),NNA4),(I(16),NNA5),(I(17),NNA6)
     1,(I(18),NNAB),(I(19),NNB1),(I(20),NNB2),(I(21),NNB3)
     1,(I(22),NNSW),(I(23),NNFGP),(I(24),NNBGP),(I(25),NNFGZ)
        EQUIVALENCE (I(26),NNBGZ),(I(27),NNB4),(I(28),NNF5),(I(29),NNB5)
     1,(I(30),NNF6),(I(31),NNB6),(I(32),NNFNP),(I(33),NNBNP)
     1,(I(34),NNBNZ),(I(35),NNBGS),(I(36),NNSR),(I(37),NNNWS)
     1,(I(38),NNG0),(I(39),NNNFC)
C     CONSTANTS NEEDED BY SUBROUTINE SETABC

C
C     SCALE FACTORS IN DB FOR GENERAL ADJUSTMENT TO:
C                    A1  A2  A3  A4  A5  A6  AN  AB  AV   AH  AF AVS
        DATA NDBSCA/-58,-65,-73,-78,-79,-80,-58,-84,-72,-102,-72,-44/
C     INCREMENT IN DB TO FORMANT AMPLITUDES OF PARALLEL BRANCH IF
C     FORMANT FREQUENCY DIFFERENCE 50, 100, 150, ... HZ
        DATA NDBCOR/10,9,8,7,6,5,4,3,2,1/
C     PRINT INPUT PAR VALUES AT T=NTIMPR, OR AT ALL TIMES IF NTIMPR=0
        DATA NTIMPR,NPPBEG,NPPEND/-1,1,39/
C
C
C   INITIALIZE SYNTHESIZER BEFORE COMPUTING WAVEFORM CHUNK IF ARG.LT.0
100     IF (INITPC.GE.0) GO TO 130
        INITPC=0
C     SET CUMULATIVE TIME COUNTER TO ZERO
        NTIMEP=0
        NAFLAS=0
C     COMPUTE SAMPLING PERIOD T (ALL CONSTANT CONTROL PARAMETERS
C     MUST BE SET BEFORE CALLING INIT)
        SAMRAT=NNSR
        T=1./SAMRAT
        PIT=PI*T
        TWOPIT=2.*PIT
        NTIMED=NNNWS*1000/NNSR
C     CONVERT INHERENTLY INTEGER PARAMS TO REAL COEFICIENTS
        C(48)=NNNWS
        C(49)=NNSW
        C(50)=NNNFC
110     CONTINUE
C
C
C
C   UPDATE ALL COEFICIENTS OF HARDWARE SYNTHESIZER
C
C     COMPUTE PARALLEL BRANCH AMPLITUDE CORRECION TO F2 DUE TO F1
130     DELF1=FLOAT(NNF1)/500.
        A2COR=DELF1*DELF1
C     COMPUTE AMPLITUDE CORRECTION TO F3-5 DUE TO F1 AND F2
        DELF2=FLOAT(NNF2)/1500.
        A2SKRT=DELF2*DELF2
        A3COR=A2COR*A2SKRT
C     TAKE INTO ACCOUNT FIRST DIFF OF GLOTTAL WAVE FOR F2
        A2COR=A2COR/DELF2
C     COMPUTE AMPLITUDE CORRECTIONS DUE TO PROXIMITY OF 2 FORMANTS
        N12COR=0
        N23COR=0
        N34COR=0
        NF21=NNF2-NNF1
        IF (NF21.LT.50) GO TO 135
        IF (NF21.LT.550) N12COR=NDBCOR(NF21/50)
        NF32=NNF3-NNF2-50
        IF (NF32.LT.50) GO TO 135
        IF (NF32.LT.550) N23COR=NDBCOR(NF32/50)
        NF43=NNF4-NNF3-150
        IF (NF43.LT.50) GO TO 135
        IF (NF43.LT.550) N34COR=NDBCOR(NF43/50)
C     PRINT INPUT PARAMETERS IF NTIMPR SET TO ZERO OR TO A SPECIFIC TIME
        IF (NTIMPR.EQ.0) GO TO 135
        IF (NTIMPR.NE.NTIMEP) GO TO 146
135     WRITE (6,140) NTIMEP
140     FORMAT (' INPUT PARS AT T = ',I4,' MS')
        WRITE (6,141) (I(NPP),NPP=NPPBEG,NPPEND)
141     FORMAT (' ',13I5)
        WRITE (6,142)
142     FORMAT (' ')
145     NPAR=1
146     NTIMEP=NTIMEP+NTIMED
C     SET AMPLITUDE OF VOICING
        NDBAV=NNG0+NNAV+NDBSCA(9)
        IMPULS=GETAMP(NDBAV)
C     AMPLITUDE OF ASPIRATION
150     NDBAH=NNG0+NNAH+NDBSCA(10)
        AHH=GETAMP(NDBAH)
C     AMPLITUDE OF FRICATION
C     (IN AN ALL-PARALLEL CONFIGURATION, AF=MAX[AF,AH])
        IF ((NNAH.GT.NNAF).AND.(NNSW.EQ.1)) NNAF=NNAH
        NDBAF=NNG0+NNAF+NDBSCA(11)
        AFF=GETAMP(NDBAF)
C     ADD A STEP TO WAVEFORM AT A PLOSIVE RELEASE
        PLSTEP=0.
        IF (NNAF-NAFLAS < 49) GO TO 151
        PLSTEP=GETAMP(NNG0+NDBSCA(11)+44_int16)
151     NAFLAS=NNAF
C     AMPLITUDE OF QUASI-SINUSOIDAL VOICING SOURCE
        NDBAVS=NNG0+NNAVS+NDBSCA(12)
        SINAMP=10.*GETAMP(NDBAVS)
C     SET AMPLITUDES OF PARALLEL FORMANTS A1 THRU A6
        NDB=NNA1+N12COR+NDBSCA(1)
        A1P=GETAMP(NDB)
        NDB=NNA2+N12COR+N12COR+N23COR+NDBSCA(2)
        A2P=A2COR*GETAMP(NDB)
        NDB=NNA3+N23COR+N23COR+N34COR+NDBSCA(3)
        A3P=A3COR*GETAMP(NDB)
        NDB=NNA4+N34COR+N34COR+NDBSCA(4)
        A4P=A3COR*GETAMP(NDB)
        NDB=NNA5+NDBSCA(5)
        A5P=A3COR*GETAMP(NDB)
        NDB=NNA6+NDBSCA(6)
        A6P=A3COR*GETAMP(NDB)
C     SET AMPLITUDE OF PARALLEL NASAL FORMANT
         NDB=NNAN+NDBSCA(7)
        ANPP=GETAMP(NDB)
C     SET AMPLITUDE OF BYPASS PATH OF FRICATION TRACT
        NDB=NNAB+NDBSCA(8)
        ABP=GETAMP(NDB)
C     RESET DIFFERENCE EQUATION CONSTANTS FOR RESONATORS
230     CALL SETABC(NNF1,NNB1,A1,B1,C1)
        CALL SETABC(NNF2,NNB2,A2,B2,C2)
        CALL SETABC(NNF3,NNB3,A3,B3,C3)
        CALL SETABC(NNF4,NNB4,A4,B4,C4)
        CALL SETABC(NNF5,NNB5,A5,B5,C5)
        CALL SETABC(NNF6,NNB6,A6,B6,C6)
        CALL SETABC(NNFNP,NNBNP,ANP,BNP,CNP)
C     AND FOR NASAL ANTIRESONATOR
        MNFNZ=-NNFNZ
        IF (MNFNZ.GE.0) MNFNZ=-1
        CALL SETABC(MNFNZ,NNBNZ,ANZ,BNZ,CNZ)
C     AND FOR GLOTTAL RESONATORS AND ANTIRESONATOR
        NPULSN=1
        IF (NNF0.LE.0) GO TO 245
C     ISSUE NO PULSE IF NNAV AND NNAVS BOTH .LE.0
        IF ((NNAV.LE.0).AND.(NNAVS.LE.0)) GO TO 245
C     WAVEFORM MORE SINUSOIDAL AT HIGH FUNDAMENTAL FREQUENCY
        NXBGP=(NNBGP*100)/NNF0
        CALL SETABC(NNFGP,NXBGP,AGP,BGP,CGP)
        CALL SETABC(0_int16,NNBGS,AGS,BGS,CGS)
        MNFGZ=-NNFGZ
        IF (MNFGZ.GE.0) MNFGZ=-1
        CALL SETABC(MNFGZ,NNBGZ,AGZ,BGZ,CGZ)
C     SET GAIN TO CONSTANT IN MID-FREQUENCY REGION FOR RGP
        AGP=.007
C     DO NOT LET F0 DROP BELOW 40 HZ
        IF (NNF0.LT.40) NNF0=40
C     MAKE AMPLITUDE OF IMPULSE INCREASE WITH INCREASING F0
        IMPULS=IMPULS*NNF0
C     NUMBER OF SAMPLES BEFORE A NEW GLOTTAL PULSE MAY BE GENERATED
        NPULSN=NNSR/NNF0
245     CONTINUE
C     CONVERT INHERENTLY INTEGER PARAMS TO REAL COEFICIENTS
        C(47)=NPULSN

        END subroutine parcoe