语音压缩

开发平台：
Unix_Linux

pctolsp2.f：源码内容
							c==========================================================================
c
c NAME	
c	pctolsp2 
c
c FUNCTION
c
c	Compute LSP from predictor polynomial.
c
c	NOTE:  Much faster conversion can be performed
c	       if the LSP quantization is incorporated.
c
c SYNOPSIS
c
c	subroutine pctolsp2(a,m,freq,lspflag)
c
c   formal 
c                       data	I/O
c	name		type	type	function
c	-------------------------------------------------------------------
c	a		real	i	a-polynomial a(0)=1
c	m		int	i	order of a
c	freq		real	o	lsp frequencies
c	lspflag		logical	o	ill-conditioned lsp test
c	n		int	na	grid points in search of zeros
c					of p-polynomials
c	eps		real	na	precision for computing zeros
c	nb		int	na	iteration limit?
c
c==========================================================================
c	
c DESCRIPTION
c
c  	Compute lsp frequencies by disection method as described in:
c	
c	Line Spectrum Pair (LSP) and Speech Data Compression,
c	F.K. Soong and B-H Juang,
c       Proc. ICASSP 84, pp. 1.10.1-1.10.4
c
c	CELP's LPC predictor coefficient convention is:
c              p+1         -(i-1)
c       A(z) = SUM   a   z          where a  = +1.0
c              i=1    i                    1
c
c	Kroon uses n=128, eps=1.e-06, nb=15 (this appears to be overkill!)
c
c**************************************************************************
c-
	subroutine pctolsp2(a,m,freq,lspflag)
c
	implicit undefined(a-z)
	integer m, maxord
	parameter (maxord=24)
	real freq(maxord), a(0:m)
	real p(maxord), q(maxord), lastfreq(maxord)
	real pi, df, ang, fm, tempfreq
	real fr, pxr, tpxr, tfr, pxm, pxl
	real fl, qxl, tqxr, tfl, qxm, qxr, tqxl
	integer mp, mh, rn, nf, mb, jc, i, j
	logical lspflag
	integer n, nb
	real eps
	parameter (n = 128)
	parameter (nb = 15)
	parameter (eps = 1.e-6)
	save lastfreq
c
	pi=atan(1.)*4.
	mp=m+1
	mh=m/2
	rn=n
c  generate p and q polynomials
	do 20 i=1,mh
	   p(i)=a(i)+a(m-i+1)
	   q(i)=a(i)-a(m-i+1)
20	continue
c  compute p at f=0.
	fl=0.
	pxl=1.
	do 30 j=1,mh
	   pxl=pxl+p(j)
30	continue
c  search for zeros of p
	nf=1
	df=.5/rn
	do 100 i=1,n
	   mb=0
	   fr=i*df
	   pxr=cos(mp*pi*fr)
	   do 80 j=1,mh
	      jc=mp-j*2
	      ang=jc*pi*fr
	      pxr=pxr+cos(ang)*p(j)
80	   continue
	   tpxr=pxr
	   tfr=fr
	   if(pxl*pxr.gt.0.)go to 95
85	   mb=mb+1
	   fm=fl+(fr-fl)/(pxl-pxr)*pxl
	   pxm=cos(mp*pi*fm)
	   do 90 j=1,mh
	      jc=mp-j*2
	      ang=jc*pi*fm
	      pxm=pxm+cos(ang)*p(j)
90	   continue
	   if(pxm*pxl.gt.0.)then
 	      pxl=pxm
	      fl=fm
	   else
	      pxr=pxm
	      fr=fm
	   endif
	   if(abs(pxm).gt.eps.and.mb.lt.4)go to 85
	   if ((pxl-pxr)*pxl .eq. 0) then
	      do 92 j=1,m
		 freq(j)=j*0.04545
92	      continue
	      print *,' pctolsp2:  default lsps used, avoiding /0'
	      return
	   end if
	   freq(nf)=fl+(fr-fl)/(pxl-pxr)*pxl
	   nf=nf+2
	   if(nf.gt.m-1)go to 120
95	   pxl=tpxr
	   fl=tfr
100	continue
c  search for the zeros of q(z)
120	continue
	freq(m+1)=.5
	fl=freq(1)
	qxl=sin(mp*pi*fl)
	do 150 j=1,mh
	   jc=mp-j*2
	   ang=jc*pi*fl
	   qxl=qxl+sin(ang)*q(j)
150	continue
	do 200 i=3,mp,2
	   mb=0
	   fr=freq(i)
	   qxr=sin(mp*pi*fr)
	   do 155 j=1,mh
	      jc=mp-j*2
	      ang=jc*pi*fr
	      qxr=qxr+sin(ang)*q(j)
155	   continue
	   tfl=fl
	   tqxl=qxl
	   tfr=fr
	   tqxr=qxr
160	   mb=mb+1
	   fm=(fl+fr)*.5
	   qxm=sin(mp*pi*fm)
	   do 180 j=1,mh
	      jc=mp-j*2
	      ang=jc*pi*fm
	      qxm=qxm+sin(ang)*q(j)
180	   continue
	   if(qxm*qxl.gt.0.)then
	      qxl=qxm
	      fl=fm
	   else
	      qxr=qxm
	      fr=fm
	   endif
	   if(abs(qxm).gt.eps*tqxl.and.mb.lt.nb)go to 160
	   if ((qxl-qxr)*qxl .eq. 0) then
	      do 93 j=1,m
		 freq(j)=lastfreq(j)
93	      continue
	      print *,' pctolsp2:  last lsps used, avoiding /0'
	      return
	   end if
	   freq(i-1)=fl+(fr-fl)/(qxl-qxr)*qxl
	   fl=tfr
	   qxl=tqxr
200	continue
c
c	*** ill-conditioned cases
c
	lspflag = .false.
	if (freq(1) .eq. 0.0 .or. freq(1) .eq. 0.5) lspflag = .true.
	do 300 i = 2, m
	   if (freq(i) .eq. 0.0 .or. freq(i) .eq. 0.5) lspflag = .true.
c			*reorder lsps if non-monotonic
	   if (freq(i) .lt. freq(i-1)) then
	      lspflag = .true.
	      print *,' pctolsp2:  non-monotonic lsps'
	      tempfreq = freq(i)
	      freq(i) = freq(i-1)
	      freq(i-1) = tempfreq
	   end if
300	continue
c			*if non-monotonic after 1st pass, reset to last values
	do 310 i=2, m
	   if (freq(i) .lt. freq(i-1)) then
	      print *,' pctolsp2:  Reset to previous lsp values'
	      do 320 j=1, m
	         freq(j) = lastfreq(j)
320	      continue
	      goto 325
	   end if
310	continue
325	do 330 i=1,m
	   lastfreq(i) = freq(i)
330	continue
	return
	end

						