home *** CD-ROM | disk | FTP | other *** search
/ Ham Radio 2000 #1 / Ham Radio 2000.iso / ham2000 / hf / dsp / source / anal.asm < prev    next >
Encoding:
Assembly Source File  |  1992-06-12  |  11.7 KB  |  499 lines
  1.     page    132,63,1,1
  2.     opt    rc
  3.     title    'LPC Analysis'
  4.  
  5. ;***************************************************************
  6. ;* ANAL.ASM -- 2400 bit/s LPC analyser module               *
  7. ;*                                   *
  8. ;* Provides pitch estimator and Leroux-Gueguen reflection      *
  9. ;* parameters estimation.                       *
  10. ;*                                   *
  11. ;* The overall implementation is based on application note     *
  12. ;*    Seshan, N.:                           *
  13. ;*    "A TMS320C30-based LPC Vocoder",               *
  14. ;*    Texas Instruments, 1990                    *
  15. ;*                                   *
  16. ;* Leroux-Gueguen algorithm implementation is based on article *
  17. ;*    LeRoux, J., Gueguen C.:                    *
  18. ;*    "A Fixed Point Computation of Partial Correlation      *
  19. ;*    Coefficients", IEEE Trans. on Acoustics, Speech and    *
  20. ;*    Signal Processing, Vol. 25, June 1977               *
  21. ;*                                   *
  22. ;* Pitch detection is based on article                   *
  23. ;*    Gold, B., Rabiner, L.:                       *
  24. ;*    "Parallel Processing Techniques for Estimating Pitch   *
  25. ;*    Periods of Speech in the Time Domain"                  *
  26. ;*    J. of the Acoustic Society of America,               *
  27. ;*    Vol. 46, August 1969                       *
  28. ;*                                   *
  29. ;* Copyright (C) 1992 by Alef Null. All rights reserved.       *
  30. ;* Author(s): Jarkko Vuori, OH2LNS                   *
  31. ;* Modification(s):                           *
  32. ;***************************************************************
  33.  
  34.     section LPCana
  35.     xdef    lpc_ana
  36.     xref    samptr
  37.     xref    p_code,g_code,k_code
  38.  
  39.     org    p:
  40.  
  41.     nolist
  42.     include 'macros'
  43.     list
  44.  
  45.  
  46. ;***************************
  47. ;*     Pitch detection       *
  48. ;***************************
  49. lpc_ana
  50. ; First lowpass filter (FIR) speech frame.
  51. ;
  52.     move            x:samptr,r0
  53.     move            #<filter+1,n0
  54.     move            #3*M-1,m0
  55.     move            #lotaps,r4
  56.     move            #<filter-1,m4
  57.     move            #temp,r1
  58.  
  59.     do    #M+2,_endlpf                ; two extra samples for pitch detector
  60.     clr    a        x:(r0)-,x0    y:(r4)+,y0
  61.     rep    #filter-1
  62.     mac    x0,y0,a     x:(r0)-,x0    y:(r4)+,y0
  63.     macr    x0,y0,a     (r0)+n0
  64.     asr    a                    ; scaling for pitch detector
  65.     move            a,x:(r1)+
  66. _endlpf
  67.  
  68. ; Generate pitch estimate for six Gold-Rabiner series
  69. ;
  70.     move            #temp,r0
  71.     move            #-1,m0
  72.     move            m0,m1
  73.     move            m1,m4
  74.     move            m4,m5
  75.     do    #M,endpeak
  76.  
  77. ; update detection envelope
  78.     move            r0,a            ; a contains speech frame index
  79.     move            #temp-1,x0
  80.     sub    x0,a        #<peakblk,r4
  81.     move            #5-1,n4
  82.     do    #6,_elope2
  83.     move            y:(r4)+,x0            ; decay if sample index > blank[j]
  84.     cmp    x0,a        y:(r4),b1
  85.     jle    <_elope1
  86.     move            y:(r4)+,x0            ; threshold[j] *= decay[j]
  87.     move            y:(r4)-,y0
  88.     mpyr    x0,y0,b
  89. _elope1 move            b1,y:(r4)+n4
  90. _elope2
  91.  
  92. ; detect peaks
  93.     move            x:(r0)+,x0
  94.     move            x:(r0)+,b
  95.     cmp    x0,b        x:(r0),x0
  96.     jle    <peak1
  97.     cmp    x0,b        (r0)-
  98.     jle    <peak2
  99.  
  100.     move            #<eblk,r1            ; update peak
  101.     move            #6,n1
  102.     move            #<peakblk,r4
  103.     move            x:(r0),b
  104.     jsr    <estimate                ; m[0]
  105.     move            x:(r0),b
  106.     move            y:<valley,y0
  107.     add    y0,b        (r1)+n1
  108.     move            (r4)+n4
  109.     jsr    <estimate                ; m[1]
  110.     move            x:(r0),b
  111.     move            y:<peak,y0
  112.     sub    y0,b        (r1)+n1
  113.     move            (r4)+n4
  114.     jsr    <estimate                ; m[2]
  115.     move            x:(r0),b
  116.     move            b,y:<peak
  117.     jmp    <peak2
  118.  
  119. ; detect valleys
  120. peak1    cmp    x0,b        (r0)-
  121.     jge    <peak2
  122.  
  123.     move            #<eblk+3*6,r1        ; update valley
  124.     move            #6,n1
  125.     move            #<peakblk+3*5,r4
  126.     move            x:(r0),b
  127.     neg    b
  128.     jsr    <estimate                ; m[3]
  129.     move            x:(r0),b
  130.     neg    b        y:<peak,y0
  131.     add    y0,b        (r1)+n1
  132.     move            (r4)+n4
  133.     jsr    <estimate                ; m[4]
  134.     move            x:(r0),b
  135.     neg    b        y:<valley,y0
  136.     sub    y0,b        (r1)+n1
  137.     move            (r4)+n4
  138.     jsr    <estimate                ; m[5]
  139.     move            x:(r0),b
  140.     neg    b
  141.     move            b,y:<valley
  142.  
  143. peak2    nop
  144. endpeak
  145.  
  146. ; update indices for next frame
  147.     move            #<peakblk,r4        ; decrement blanking indices
  148.     move            #5,n4
  149.     move            #<peakblk+3,r5        ; decrement series indices
  150.     move            n4,n5
  151.     move            #>M,x0
  152.     do    #6,_updpk
  153.     move            y:(r4),a
  154.     sub    x0,a        y:(r5),b
  155.     sub    x0,b        a,y:(r4)+n4
  156.     move            b,y:(r5)+n5
  157. _updpk
  158.  
  159. ; Choose best pitch from Gold-Rabiner 6*6 pitch estimate matrix
  160. ;
  161.     clr    a
  162.     move            a,y:<p_a
  163.     move            #<voiced,a1         ; assume frame is unvoiced initially
  164.     move            a1,y:<c_best
  165.  
  166. ; calculate last columns of pitch estimate matrix
  167.     move            #<eblk,r1
  168.     move            #2,n1
  169.     do    #6,_cpch2
  170.     do    #3,_cpch1
  171.     move            y:(r1)+,a            ; e[i][j+3] = e[i][j]+e[i][j+1]
  172.     move            y:(r1)+,x0
  173.     add    x0,a        (r1)+
  174.     move            a,y:(r1)-n1
  175. _cpch1    move            (r1)+n1
  176.     move            (r1)+
  177. _cpch2
  178.  
  179. ; use first column of pitch estimates as potential cadidates
  180.     move            #<eblk,r1
  181.     move            #6,n1
  182.     do    #6,pold1
  183.  
  184. ; compute coincidence count under four tolerance levels
  185.     move            #<bias,r4
  186.     do    #4,pold2
  187.     move            y:(r4)+,b            ; c = bias[k]
  188.     move            b,y:<c
  189.  
  190.     move            #>bias,x0            ; a = 0.04*e[l][0]*(k+1)+.5
  191.     move            r4,b
  192.     sub    x0,b        y:(r1),x0
  193.     move            b,x1
  194.     mpy    x0,x1,b     #0.04,x0            ; integer to integer multiplication
  195.     asr    b                    ; adjust binary point
  196.     move            b0,x1
  197.     mpyr    x0,x1,a     #<eblk,r0            ; integer to fraction multiplication
  198.  
  199. ; compare against all thirty-six values in estimate matrix
  200.     do    #6*6,_pcomp2
  201.  
  202. ; if compare estimate is within tolerance increment coincidence count
  203.     move            y:(r1),b            ; if abs(e[l][0]-e[i][j]) < a
  204.     move            y:(r0)+,y0
  205.     sub    y0,b
  206.     abs    b        #>1,y0
  207.     cmp    b,a        y:<c,b
  208.     jlt    <_pcomp1
  209.     add    y0,b                    ; then increment coincidence count
  210.     move            b,y:<c
  211. _pcomp1 nop
  212. _pcomp2
  213.  
  214. ; if coincidence count is higher then for previous best estimate current
  215. ; candidate is new best estimate and current coincidence count is best coincidence count
  216.     move            y:<c_best,a
  217.     cmp    a,b        y:(r1),a
  218.     jlt    <pold
  219.     move            b,y:<c_best
  220.     move            a,y:<p_a
  221. pold    nop
  222. pold2    move            (r1)+n1
  223. pold1
  224.  
  225. ; Post process pitch
  226. ; Correct possible pitch tracking/voicing errors
  227. ;
  228.  
  229. ; if an unvoiced frame occurs between two voiced frames interpolate
  230.     move            y:<p_a,a            ; if p_a && !p_a_0 && p_a_1
  231.     tst    a        y:<p_a_0,b
  232.     jeq    <voice
  233.     tst    b        y:<p_a_1,a
  234.     jne    <voice
  235.     tst    a        y:<p_a,b
  236.     jeq    <voice
  237.  
  238.     add    b,a                    ; then p_code = (p_a + p_a_1)/2
  239.     asr    a
  240.     rnd    a
  241.     jmp    <median
  242.  
  243. ; else if a voiced frame occurs between two unvoice frames switch to unvoice
  244. voice    move            y:<p_a,a
  245.     tst    a        y:<p_a_0,b
  246.     jne    <nope
  247.     tst    b        y:<p_a_1,a
  248.     jeq    <nope
  249.     tst    a
  250.     jne    <nope
  251.  
  252.     clr    a
  253.     jmp    <median
  254.  
  255. nope    move            y:<p_a_0,a
  256.  
  257. ; perform median filtering within a series of voiced frames
  258. median    move            a,x:<p_code
  259.  
  260.     move            y:<p_a,a
  261.     tst    a        y:<p_a_0,b
  262.     jeq    <medend
  263.     tst    b        y:<p_a_1,a
  264.     jeq    <medend
  265.     tst    a        x:<p_code,b
  266.     jeq    <medend
  267.  
  268.     move            y:<p_a,a            ; i=a,j=b,l=x0,k=y0
  269.     move            b,x0
  270.  
  271.     cmp    x0,a        y:<p_a_1,y0         ; j = (i>l) ? i : j
  272.     tgt    a,b
  273.     cmp    y0,b                    ; j = (j>k) ? k : j
  274.     tgt    y0,b
  275.     cmp    x0,a                    ; i = (i>l) ? l : i
  276.     tgt    x0,a
  277.     cmp    b,a                    ; b = (i>j) ? i : j
  278.     tgt    a,b
  279.     move            b,x:<p_code
  280.  
  281. medend    move            y:<p_a_0,a1         ; update previous pitch locations
  282.     move            a1,y:<p_a_1
  283.  
  284.     move            y:<p_a,a1
  285.     move            a1,y:<p_a_0
  286.  
  287. ;***************************
  288. ;* Preemphasize and window *
  289. ;***************************
  290.  
  291. ; Preemphasize (H(z) = 1 - az^-1) speech signal and apply hamming window
  292.     move            x:samptr,r0
  293.     move            #temp,r1
  294.     move            #hamming,r4
  295.     move            #-0.9375,x1         ; a = 15/16 = 0.9375
  296.  
  297.     move            x:(r4)+,y0            ; first term
  298.     move            x:(r0)+,x0
  299.     mpyr    x0,y0,a
  300.     move            a,x:(r1)+
  301.  
  302.     do    #N-1,_hwin                ; remaining terms
  303.     move            x:(r0),a
  304.     macr    x0,x1,a     x:(r4)+,y0
  305.     move            a,y1
  306.     mpyr    y0,y1,a     x:(r0)+,x0
  307.     move            a,x:(r1)+
  308. _hwin
  309.  
  310.  
  311. ;***************************
  312. ;*     Autocorrelation       *
  313. ;***************************
  314.  
  315. ; Returns short time autocorrelation of result
  316.     move            #temp,r0
  317.     move            #en,r4
  318.     move            #<N,r5
  319.  
  320.     do    #P+1,_endac
  321.     move            #temp,r2
  322.     move            r0,r1
  323.     clr    a        x:(r2)+,x1
  324.     do    r5,_endlag
  325.     move            x:(r1)+,x0
  326.     mac    x0,x1,a     x:(r2)+,x1
  327. _endlag
  328.     rnd    a        (r5)-
  329.     move            x:(r0)+,a    a,y:(r4)+   ; just incrementing r0
  330. _endac
  331.  
  332. ; Setup ACs for LG-algorithm
  333.     move            #en+1,r0
  334.     move            #ep,r1
  335.  
  336.     do    #P,_accpy
  337.     move            y:(r0)+,a1
  338.     move            a1,y:(r1)+
  339. _accpy
  340.  
  341.  
  342. ; store previous results (in order to compensate
  343. ; pitch post processing delay)
  344.     move            y:g_a,a1
  345.     move            a1,x:g_code
  346.  
  347.     move            #k_a,r0
  348.     move            #k_code,r1
  349.     do    #P,_endsto
  350.     move            y:(r0)+,a1
  351.     move            a1,x:(r1)+
  352. _endsto
  353.  
  354. ;***************************
  355. ;*    Leroux-Gueguen       *
  356. ;***************************
  357.  
  358. ; Derive reflection coefficients using Leroux-Gueguen algorithm
  359.     move            #ep,r0
  360.     move            #k_a,r4
  361.     do    #P,lgloop
  362.  
  363.     move            y:en,x0            ; k_a[j] = -ep[j]/en[0]
  364.     move            y:(r0),a
  365.     abs    a        a,b
  366.     eor    x0,b        b,y0
  367.     and    #$fe,ccr
  368.     rep    #24
  369.     div    x0,a
  370.     jmi    _L1
  371.     neg    a
  372. _L1    move            a0,y:(r4)
  373.  
  374.     move            r0,r5
  375.     move            #en,r1
  376.     move            y:(r4)+,x0
  377.     do    lc,_lgker
  378.  
  379.     move            y:(r5),x1            ; temp     = en[i]
  380.     move            y:(r1),a            ; en[i]   += k_a[j]*ep[i+j]
  381.     macr    x0,x1,a     y:(r1),y0            ; ep[i+j] += k_a[k]*temp
  382.     move            x1,b
  383.     macr    x0,y0,b     a,y:(r1)+
  384.     move            b,y:(r5)+
  385. _lgker
  386.     move            (r0)+
  387. lgloop
  388.     move            y:en,x0            ; g_a = sqrt(en[0]/N)
  389.     move            #1.0/N,x1
  390.     mpyr    x0,x1,a
  391.     move            a,y:g_a
  392.  
  393.     rts
  394.  
  395.  
  396. ; Estimate peaks by checking if potential series value is outside blanking
  397. ; interval and exceeds threshold value. If a new pitch estimate is generated
  398. ; and Gold-Rabiner parameters are updated.
  399. ;   r1 - pointer to e vector (a row of e matrix)
  400. ;   r4 - pointer to (blank,threshold,decay,i0,p_av) vector (incremented by 1)
  401. ;    a - speech frame sample index
  402. ;    b - Gold-Rabiner series value
  403. estimate move            y:(r4)+,x0            ; if index > blank[j] and value > threshold[j]
  404.     cmp    x0,a        y:(r4),x0
  405.     jle    <estiend
  406.     cmp    x0,b
  407.     jle    <estiend
  408.  
  409. ; update envelope threshold
  410.     tfr    a,b        b,y:(r4)+
  411.  
  412. ; update 3rd and 2nd newest estimates
  413.     move            y:(r1)+,x0            ; e[j][2] = e[j][1]
  414.     move            y:(r1)+,x1            ; e[j][1] = e[j][0]
  415.     move            x1,y:(r1)-
  416.     move            x0,y:(r1)-
  417.  
  418. ; compute newest estimate and update index to last point of pitch estimation
  419.     move            (r4)+
  420.     move            y:(r4),x0           ; e[j][0] = index - i0[j]
  421.     move            b,y:(r4)+           ; i0[j] = index
  422.     sub    x0,b        y:(r4),x0
  423.     move            b,y:(r1)
  424.  
  425. ; compute moving pitch average
  426.     add    x0,b        #>32,x0            ; p_av[j] = 0.5*(e[j][0]+p_av[j])
  427.     asr    b        #>80,x1
  428.     rnd    b        #<exptab,r2
  429.     cmp    x0,b                    ; limit to interval 32..80
  430.     tlt    x0,b
  431.     cmp    x1,b        #0.4,y1
  432.     tgt    x1,b
  433.     sub    x0,b        b,y:(r4)
  434.  
  435. ; calculate envelope decay value
  436.     move            b,n2
  437.     move            y:(r4)-,x1
  438.     move            y:(r2+n2),y0        ; decay[j] = exp(-.695/p_av[j])
  439.  
  440. ; update blanking interval
  441.     mpy    x1,y1,b     (r4)-            ; blank[j] = .4*p_av[j]+i
  442.     add    a,b        y0,y:(r4)-
  443.     rnd    b        (r4)-
  444.     move            b,y:(r4)+
  445.  
  446. estiend rts
  447.  
  448.  
  449. ;****************************
  450. ;*     DATA - AREA        *
  451. ;****************************
  452.  
  453.     org    y:
  454.  
  455. ; Gold-Rabiner pitch tracking variables
  456. p_a    ds    1                    ; pitch estimate for current frame
  457. p_a_0    dc    0                    ; previous frame's p_a
  458. p_a_1    dc    0                    ; previous frame's p_a_0
  459. c_best    ds    1
  460.  
  461. c    ds    1                    ; coincidence counter
  462.  
  463. peak    dc    0                    ; last maximum in speech frame
  464. valley    dc    0                    ; last minimum in speech frame
  465.  
  466. peakblk dc    0,0,0,0,32                ; blank, threshold, decay, i0, p_av
  467.     dc    0,0,0,0,32
  468.     dc    0,0,0,0,32
  469.     dc    0,0,0,0,32
  470.     dc    0,0,0,0,32
  471.     dc    0,0,0,0,32
  472.  
  473. eblk    dc    0,0,0,0,0,0                ; pitch estimate matrix
  474.     dc    0,0,0,0,0,0
  475.     dc    0,0,0,0,0,0
  476.     dc    0,0,0,0,0,0
  477.     dc    0,0,0,0,0,0
  478.     dc    0,0,0,0,0,0
  479.  
  480. bias    dc    -2,-3,-6,-8                ; tolerance based coincidence bias
  481. exptab                            ; decay exponent table
  482. dum    set    32
  483.     dup    80-32+1
  484.     dc    @xpn(-.695/@cvf(dum))
  485. dum    set    dum+1
  486.     endm
  487.  
  488. ; Leroux-Gueguen coefficient extraction variables
  489. k_a    ds    P                    ; reflection coefficient array
  490.  
  491. en    ds    P+1                    ; Leroux-Gueguen estimation arrays
  492. ep    ds    P
  493.  
  494. g_a    dc    0                    ; gain
  495.  
  496.     endsec
  497.  
  498.     end
  499.