开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 通讯/手机编程 > 语音压缩 > 查看源码
celp.F
资源名称：celp_3.2a.tar.gz　[点击查看]
上传用户：szhypcb168
上传日期：2007-01-06
资源大小：2187k
文件大小：32k
源码类别：

语音压缩

开发平台：

Unix_Linux

celp.F：源码内容
  1. ***************************************************************************
  2. *                                                                         *
  3. * CELP Voice Coder                                                  *
  4. * Version 3.2                                                   *
  5. *                                                                         *
  6. * The U.S. Government shall not be held liable for any damages      *
  7. * resulting from this code.  Further reproduction or distribution   *
  8. * of this code without prior written permission of the U.S.         *
  9. * Government is prohibited.  See NCS Technical Information Bulletin *
  10. * on FS 1016 CELP Speech Coding.   *                                                                         *   *
  11. ***************************************************************************
  12. C
  13. C ROUTINE
  14. C celp (main)
  15. C
  16. C FUNCTION
  17. C Code excited linear predictor (main routine)
  18. C
  19. C==========================================================================
  20. C
  21. C REFERENCES
  22. C
  23. C National Communication System Technical Information Bulletin
  24. C Federal Standard 1016 (to be published 1992).
  25. C
  26. C Campbell, Joseph P. Jr., Thomas E. Tremain and Vanoy C. Welch,
  27. C "The Federal Standard 1016 4800 bps CELP Voice Coder," Digital
  28. C Signal Processing, Academic Press, Vol1, No. 3, p. 145-155.
  29. C
  30. C Kemp, David, P., Retha A. Sueda and Thomas E. Tremain, "An
  31. C Evaluation of 4800 bps Voice Coders," Proceedings of the IEEE
  32. C International Conference on Acoustics, Speech and Signal Processing
  33. C (ICASSP), 1989, p. 200-203.
  34. C
  35. C Fenichel, R., "Federal Standard 1016," National Communications
  36. C System, Office of Technology and Standards, Washington, DC 20305-2010,
  37. C 14 February 1991.
  38. C
  39. C Campbell, Joseph P. Jr., Thomas E. Tremain and Vanoy C. Welch,
  40. C "The DoD 4.8 kbps Standard (Proposed Federal Standard 1016),"
  41. C "Advances in Speech Coding", Kluwer Academic Publishers, 1991,
  42. C Chapter 12, p. 121-133.
  43. C
  44. C   Tutorials:
  45. C Fallside, Frank and William Woods, Computer Speech Processing,
  46. C Prentice Hall International, 1985, Chapter 4 (by Bishnu Atal).
  47. C
  48. C==========================================================================
  49. C
  50. C 4800 bps CELP Characteristics
  51. C
  52. C                  Spectrum       Pitch            Code Book
  53. C                  -------------  ---------------  ----------------
  54. C       Update     30 ms          30/4 = 7.5 ms    30/4 = 7.5 ms
  55. C                  ll=240         lp=60            l=60
  56. C
  57. C       Order      10             256 (max) x 60   512 (max) x 60
  58. C                                 1 gain           1 gain
  59. C
  60. C       Analysis   Open loop      Closed loop      Closed loop
  61. C                  Correlation    Modified MSPE    MSPE VQ
  62. C                  30 ms Hamming  VQ, weight=0.8   weight=0.8
  63. C                  no preemphasis range 20 to 147  shift by 2
  64. C                  15 Hz BW exp   (w/ fractions)   77% sparsity
  65. C
  66. C       Bits per   34 indep LSP   index:  8+6+8+6  index:     9*4
  67. C        Frame     [3444433333]   gain(-1,2): 5*4  gain(+/-): 5*4
  68. C
  69. C       Bit Rate   1133.3 bps     1600 bps         1866.67 bps
  70. C
  71. C NOTE:  The remaining 200 bps are used as follows:  1 bit per frame
  72. C for synchronization, 4 bits per frame for forward error correction
  73. C and 1 bit per frame to provide future expansion(s) of the coder.
  74. C
  75. C==========================================================================
  76. C
  77. C UNPERMUTED BIT ASSIGNMENT
  78. C
  79. C lsp  1 1-3 lsp  6 20-22
  80. C lsp  2 4-7 lsp  7 23-25
  81. C lsp  3 8-11 lsp  8 26-28
  82. C lsp  4 12-15 lsp  9 29-31
  83. C lsp  5 16-19 lsp 10 32-34
  84. C
  85. C Subframe:   1   2   3    4
  86. C ---------- ----- ----- -----  -------
  87. C       pitch delay 35-42 ..... 87-94  .......
  88. C       delta delay ..... 62-67 .......  114-119
  89. C       pitch gain  43-47 68-72 95-99  120-124
  90. C       cbindex 48-56 73-81 100-108  125-133
  91. C       cbgain 57-61 82-86 109-113  134-138
  92. C
  93. C       future bit 139
  94. C       error control 140-143
  95. C       sync 144
  96. C
  97. C FORWARD ERROR CONTROL, HAMMING (15,11):
  98. C bits protected allocation rate parameter
  99. C 0 1 LSP
  100. C 40,41,42,92,93,94 1st-3rd MSBs 2 Pitch delay
  101. C 0 2 Delta pitch delay
  102. C 47,72,99,124     1st MSB 4 Pitch gain
  103. C 0 4 CB index
  104. C 0 4 CB gain
  105. C 139 1 1 Future bit
  106. C
  107. C SYNCHRONIZATION BIT
  108. C The sync bit (144) begins with 0 in the first frame, then alternates 
  109. C between 1 and 0 on successive frames.
  110. C
  111. C==========================================================================
  112. C 
  113. C INPUT FILES
  114. C ifile.spd  speech
  115. C ifile.chan channel (use for synthesis only)
  116. C lsp34.table spectrum (34 bit LSP) coding
  117. C pgain.table pitch gain gain coding
  118. C cbgain.table code book gain coding
  119. C codebook.h stochastic code book
  120. C pdelay.h adaptive code book (pitch) delay
  121. C pdencode.h pitch delay bit assignment encoding table
  122. C pddecode.h pitch delay bit assignment decoding table
  123. C submult.h submultiple pitch delay pointer table
  124. C bitprot.h bit protection (specifies FEC bits)
  125. C bitperm.h bit permutation (for burst errors)
  126. C
  127. C OUTPUT FILES
  128. C ofile.spd  speech (without high pass)
  129. C ofilehpf.spd speech (high passed)
  130. C ofilenpf.spd speech (no post filtering)
  131. C ofile.chan bit stream channel file
  132. C
  133. #ifdef SUNGRAPH
  134. C Sungraph (tm) Files:
  135. C ifile_ofile.sg_data  
  136. C channel.sg_data  
  137. C stream_error.sg_data
  138. C lsp1.sg_data  see sungraph_open.com
  139. C lsp2.sg_data
  140. C codebook.sg_data
  141. C constrain.sg_data
  142. C pitch.sg_data
  143. C error.sg_data
  144. C rc.sg_data
  145. #endif
  146. C**************************************************************************
  147. C*-
  148.         program celp
  149. implicit undefined (a-z)
  150. include 'ccsub.com'
  151. convex #include "ccsub.com"
  152.         real sold(maxll), snew(maxll), ssub(maxll)
  153. real v(maxll), vdecoded(maxll), hamw(maxll), hamws(maxll)
  154. real dpa(maxpa), dps(maxpa), ahpf(3), bhpf(3)
  155. real dppa(maxpa), dpps(maxpa)
  156. real ahpfo(3), bhpfo(3)
  157. real newfreq(maxno), unqfreq(maxno), rcn(maxno)
  158. real lsp(maxno,maxll/maxl), scale, descale, decodedgain, taus(4)
  159. real lambda, omega, alpha, beta, pddecode(0:maxpd-1)
  160. real predgain, preddiff, predavg, mindiff, dm(9), dm2(9)
  161. real cbg(maxll/maxl), pgs(maxll/maxl)
  162. integer cbi(maxll/maxl), pdencode(0:maxpd-1), pdtabi(0:maxpd-1)
  163. integer framesnr, framedm, framedm2, lnblnk
  164. integer i, j, k, l, ll, lp, nn, np
  165. integer nrec1, nrec2, nrec3, nrec9,  ierr, iodisk, oflen
  166. integer findex(maxno), error, total, ranseed, iflen
  167. integer*2 iarf(maxll), npf(maxll), pf(maxll)
  168. logical flag, snrflag, lspflag, analy, clip
  169. character*80 ifile, ofile, chfile, mfile
  170. character stype*10, anal*10
  171. data ranseed /1234567/, frame /1/, framesnr/0/
  172. data nrec1 /1/, nrec2 /1/, nrec3 /1/, nrec9 /1/
  173. c
  174. c *bit stream
  175. integer cbbits, pointer, streambits, ssum, mask(172)
  176. integer bitsum1, bitsum2, bitpointer, sbits(maxno), sync
  177. data sync /1/
  178. parameter (streambits=144)
  179. integer*2 stream(streambits), oldstream(streambits)
  180. integer pstream(streambits)
  181. integer*2 savestream(streambits), iint
  182. character*36 line
  183. c
  184. c *filter memories
  185. integer maxnop1
  186. parameter (maxnop1=maxno+1)
  187. real dhpf1(3), dhpf2(3), dsa(maxno+1), dss(maxno+1)
  188. real dhpf1o(3), dhpf2o(3)
  189. real dp1(maxno+1), dp2(maxno+1), dp3(2), dp4(2)
  190. real ip, op, sumsnr, sumdm(10), sumdm2(10)
  191. c
  192. c *initialize filter memories
  193. data d2a/maxnop1*0./, d2b/maxnop1*0./ 
  194. data d3a/maxnop1*0./, d3b/maxnop1*0./
  195. data d4a/maxnop1*0./, d4b/maxnop1*0./, d5/maxnop1*0./
  196. data d1a/maxpa*0./, d1b/maxpa*0./, d6/maxpa*0./
  197. data dhpf1/3*0./, dhpf2/3*0./, dsa/maxnop1*0./, dss/maxnop1*0./
  198. data dhpf1o/3*0./, dhpf2o/3*0./, syndavg/0./
  199. data dp1/maxnop1*0./, dp2/maxnop1*0./, dp3/2*0./, dp4/2*0./
  200. data ip/0./, op/0./, sumsnr/0./, sumdm/10*0/, sumdm2/10*0/
  201. data dpa/maxpa*0./, dps/maxpa*0./,dppa/maxpa*0./, dpps/maxpa*0./
  202. c
  203. c *error control coding parameters:
  204. real syndrun, ber, realerror, syndavg
  205. parameter (syndrun = 100.)
  206. integer codelength1, codelength2, paritylength
  207. parameter (codelength1 = 15)
  208. parameter (codelength2 = 11)
  209. parameter (paritylength = codelength1 - codelength2)
  210.         integer bitprotect(codelength2), codeword(codelength1)
  211.         integer hmatrix(codelength1), syndrometable(codelength1)
  212. integer paritybit, errorflag, bitpermute(streambits)
  213. integer perrorcount, syndrome, eccbits, errorcount
  214. logical twoerror, protect
  215. c
  216. #ifdef SUNGRAPH
  217. c
  218. c *** sungraph
  219. c *define sungraph functions
  220.         integer end_block, close_file, status
  221.         integer open_file, def_variable
  222. integer close_channel, open_var_channel, read_variable
  223. integer input_fid, num_read
  224. c *define sungraph files
  225. integer ifile_ofile_fid, channel_fid, lsp1_fid, lsp2_fid, constrain_fid
  226. integer cb_fid, pitch_fid, error_fid, stream_error_fid, rc_fid
  227. character temp_fname*80, str*17
  228. c *symbolic disk_io data types
  229. integer i2, i4, r4
  230. parameter (i2 = 1)
  231. parameter (i4 = 2)
  232. parameter (r4 = 3)
  233. c *nonhigh-passed old input speech for alignment
  234. real ssubnhp(maxll), soldnhp(maxll)
  235. c *sungraph variable id's
  236. include 'sgvar.com'
  237. convex #include "sgvar.com"
  238. #endif
  239. c
  240. c *** defaults
  241. c
  242. c Default command line i/o file names:
  243. data ifile /'../speech/3m3f.spd'/
  244. data ofile /'ofile'/, chfile /'chfile'/
  245. c Default code book parameters:
  246. data ncsize /512/, cbbits /9/, l /60/, mxsw /.true./
  247. c Default LPC analysis parameters:
  248. data ll /240/, no /10/
  249. c Default pitch analysis parameters:
  250. data lp /60/, np /1/, pstype /'hier'/, prewt /0.0/
  251. c Default noise weighting factor:
  252. data gamma /0.8/
  253. c Default input speech scale factor:
  254. data scale /1.0/, descale /1.0/
  255. c Default code book gain quantization:
  256. data cbgtype /'log'/, cbgbits /5/
  257. c Default pitch bit allocation (tau, np a's) (maxnp+1-2=2):
  258. data ptype /'max2'/, pbits /8, 6, 5, 2*0/
  259. c Default spectrum bit allocation (maxno-10=14):
  260. data stype /'kang'/, sbits /3,4,4,4,4,3,3,3,3,3/
  261. c Default error rate (eccbits=fecbits+future bit)
  262. data ber /0.0/, error /0/, total /0/, eccbits /5/
  263. c Default mask file
  264. data mfile /'none'/
  265. c
  266. c Spectrum analysis parameters (should be in cli!)
  267. data anal /'autohf'/
  268. copt Filter Coefficients for 2nd order 100 Hz HPF with 60 Hz notch:
  269. copt filter data ahpf/1., -1.99778, 1./, bhpf/1., -1.88, .89/
  270. c Filter Coefficients for 2nd order Butterworth 100 Hz HPF:
  271. data ahpf/0.946, -1.892, 0.946/, bhpf/1., -1.889033, 0.8948743/
  272. c Filter Coefficients for 2nd order Butterworth 275 Hz HPF:
  273. c275 data ahpfo/0.858, -1.716, 0.858/, bhpfo/1., -1.696452, 0.7368054/
  274. c Filter Coefficients for 2nd order Butterworth 100 Hz HPF:
  275. data ahpfo/0.946, -1.892, 0.946/, bhpfo/1., -1.889033, 0.8948743/
  276. c High frequency (stabilization) factor
  277. parameter (lambda = 0.0)
  278. c Bandwidth expansion for LPC analysis (15 Hz)
  279. parameter (omega = 0.994127)
  280. c Bandwidth expansions for postfilter
  281. parameter (alpha = 0.8)
  282. parameter (beta  = 0.5)
  283. c
  284. c *** handle underflows on Sun (Sun3 or Sun4)
  285. c This issue can be reviewed by reading the UNIX man pages
  286. c on ieee floating point (man -k ieee).
  287. c
  288. call abrupt_underflow()
  289. #include <f77/f77_floatingpoint.h>
  290. cSUN? call ieee_handler("set", "underflow", SIGFPE_DEFAULT)
  291. c
  292. c *** parse the command line
  293. c
  294. c celp [-i ifile] [-o ofile] [-c chan]
  295. c      [-p pfile] [-q qfile] [-m mfile] [-l lfile]
  296. c
  297. c note:  ncsize, no, and gamma are passed through common
  298. c note:  The calls to cli & cliend may be commented out and celp will
  299. c        operate properly with default values (without bells & whistles)
  300. c
  301. c ... If synthesizer only, 
  302. c
  303. c ... If synthesizer only, the command line interface (cli) is not
  304. c     used.  (Separate analysis is not provided since the combined
  305. c     analyzer and synthesizer generates a bit stream channel file.)
  306. c     The celp command syntax for synthesis only (analy=.false.) is:
  307. c
  308. c celp [input.chan] [ofile]
  309. c
  310. c where: input is input.chan ascii hex bit stream channel file
  311. c        output is ofile.spd speech file (postfiltered)
  312. c
  313. c *read command line,
  314. analy = .true.
  315. call cli(ifile, ofile, l, ll, lp, np, scale, descale, 
  316.      &           ber, mask, stype, sbits, eccbits, ssum, analy)
  317. c
  318. c *** initialize
  319. c *number of code words/LPC frame
  320. nn = ll/l
  321. c *dimension of d1a, d1b and d6???
  322. idb = mmax+maxnp-1+l
  323. c
  324. plevel1 = 2**pbits(1)
  325. c *levels of delta tau
  326. plevel2 = 2**pbits(2)
  327. c
  328. c *number of bits per subframe
  329. bitsum1 = cbbits+cbgbits+pbits(1)+pbits(3)
  330. bitsum2 = cbbits+cbgbits+pbits(2)+pbits(3)
  331. c
  332. c *enable/disable error control coding
  333. protect  = .true.
  334. errorflag= 0
  335. c *for double error detecting FEC codes (NOT USED)
  336. twoerror = .false.
  337. c
  338. snrflag  = .false.
  339. lspflag  = .true.
  340. c
  341. c *start nseg at 0 to do pitch on odd segments
  342. c (nseg is incremented before csub) 
  343.         nseg  = 0
  344. c
  345. c *generate matrix for error control coding
  346. call matrixgen(codelength1, codelength2, hmatrix, syndrometable)
  347. c
  348. c *generate Hamming windows
  349.         call ham(hamw,ll)
  350. c *UNNECESSARY, used for distortion diagnostics
  351.         call ham(hamws,l)
  352. c
  353. c
  354. c *** open and define files
  355. c
  356. if (analy) then
  357. #ifdef SUNGRAPH
  358. c
  359. c *sungraph (tm) files
  360. include 'sgopen.com'
  361. convex #include "sgopen.com"
  362. #else
  363. c *input file
  364. if (iodisk(3, 9, ifile, nrec9, iarf, ll) .ne. 0)
  365.      &    stop '*** Error opening ifile.spd'
  366. #endif
  367. end if
  368. c *postfiltered & nonpostfiltered output
  369. oflen = lnblnk(ofile)
  370. iflen = lnblnk(ifile)
  371. if(iodisk(4,1,ofile(1:oflen)//'.spd',nrec1,pf,l).ne.0)
  372.      &    stop '*** Error opening ofile.spd'
  373. cnpf if(iodisk(4,2,ofile(1:oflen)//'npf.spd',nrec2,npf,l).ne.0)
  374. cnpf     +    stop '*** Error opening ofilenpf.spd'
  375. if(iodisk(4,3,ofile(1:oflen)//'hpf.spd',nrec3,pf,l).ne.0)
  376.      &    stop '*** Error opening ofilehpf.spd'
  377. c
  378. c *bit stream channel file
  379. if (analy) then
  380.    open(unit=25, file=ofile(1:oflen)//'.chan')
  381. else
  382.    open(unit=25, file=ifile(1:iflen)//'.chan', status='old')
  383. end if
  384. c
  385. c *quantizer design (data collection) files
  386. copt open(unit=20, file='cbgain.data')
  387. copt open(unit=21, file='pgain.data')
  388. c
  389. c *read stochastic code book vector file (by 20)
  390. open(unit=11, file='./codebook.h', status='old', err=112)
  391. i = 1
  392. do 110 j = 1, (2*(maxncsize-1)+l)/20
  393.    read(11, *, err=112) (x(k), k=i,i-1+20)
  394.    i = i + 20
  395. 110 continue
  396. read(11, *, err=112) (x(k), k=i,i-1+mod((2*(maxncsize-1)+l),20))
  397. close(11)
  398. go to 115
  399. 112    stop ' celp:  Problem with file "codebook.h"'
  400. 115 continue
  401. c
  402. c *read adaptive code book index (pitch delay) file
  403. open(unit=12, file='./pdelay.h', status='old', err=122)
  404. do 120 i = 0, maxpd-1
  405.    read(12, *, err=122) pdelay(i)
  406. 120 continue
  407. close(12)
  408. go to 125
  409. 122    stop ' celp:  Problem with file "pdelay.h"'
  410. 125 continue
  411. c
  412. c *read pitch delay coding tables for bit assignment
  413. c *pdencode.h for encoding, pddecode.h for decoding
  414. c *generate pdtabi for delta delay coding
  415. open(unit=13, file='./pdencode.h', status='old', err=132)
  416. open(unit=14, file='./pddecode.h', status='old', err=142)
  417. do 130 i = 0, maxpd-1
  418.    read(13, '(z)', err=132) pdencode(i)
  419.    read(14, *, err=142) pddecode(i)
  420.    pdtabi(pdencode(i)) = i
  421. 130 continue
  422. close(13)
  423. close(14)
  424. go to 145
  425. 132    stop ' celp:  Problem with file "pdencode.h"'
  426. 142    stop ' celp:  Problem with file "pddecode.h"'
  427. 145 continue
  428. c
  429. c *read bit protection file
  430. open(unit=15, file='./bitprot.h', status='old', err=152)
  431.    read(15, *, err=152) bitprotect(1:codelength2)
  432. close(15)
  433. go to 155
  434. 152    stop ' celp:  Problem with file "bitprot.h"'
  435. 155 continue
  436. c
  437. c *read bit permutation file (by 12)
  438. open(unit=16, file='./bitperm.h', status='old', err=162)
  439. i = 1
  440. do 160 j = 1, streambits/12
  441.    read(16, *, err=162) (bitpermute(k), k=i,i+11)
  442.    i = i + 12
  443. 160 continue
  444. close(16)
  445. go to 165
  446. 162    stop ' celp:  Problem with file "bitperm.h"'
  447. 165 continue
  448. c
  449. c..................... m a i n  l o o p ........(fortran 77 do while)......
  450. c
  451. c
  452. c *** ANALYSIS ..................................................
  453. c
  454. c *** if synthesizer only, skip analyzer >>>>>>>>>>>>>>>>>>>>>>>>
  455. c
  456. 69 if (analy) then
  457.    pointer = 1
  458. c
  459. c *read speech segment s of size ll, until end of file
  460. c
  461. #ifdef SUNGRAPH
  462.    status=read_variable(input_fid, iarf, ll, num_read)
  463.    if (status .lt. 0) then
  464.       if (status .eq. -5) go to 999
  465.       call print_disk_read_error(status,'read_input_fid')
  466.    end if
  467. #else
  468.    if (iodisk(1, 9, ifile, nrec9, iarf, ll) .ne. ll) go to 999
  469. #endif
  470. c
  471. c *scale and convert to real speech 
  472. c *The ssub buffer used for subframe CELP analysis is 1/2 a
  473. c *frame behind the snew buffer and 1/2 a frame ahead of the 
  474. c *sold buffer.
  475.    do 200 i=1,ll
  476.               snew(i)=max(min(scale * iarf(i), 32767.), -32768.)
  477. 200    continue
  478. #ifdef SUNGRAPH
  479. c
  480. c *create ssubnhp vector for sungraph
  481.    do 203 i=1,ll/2
  482.       ssubnhp(i)=soldnhp(i+ll/2)
  483.       ssubnhp(i+ll/2)=snew(i)
  484. 203    continue
  485. c
  486. c *save input speech in file 'ifile_ofile'
  487.    call save_sg(ifile_vid, ssubnhp, ll, 'save speech_in_vid')
  488. c
  489. c *save snew in soldnhp for sungraph in next frame
  490.    do 207 i=1,ll
  491.       soldnhp(i)=snew(i)
  492. 207    continue
  493. #endif
  494. c
  495. c *high pass filter snew
  496.    call zerofilt(ahpf,2,dhpf1,snew,ll)
  497.    call polefilt(bhpf,2,dhpf2,snew,ll)
  498. c
  499. c
  500. c *make ssub vector from snew and sold
  501.    do 210 i=1,ll/2
  502.       ssub(i)=sold(i+ll/2)
  503.       ssub(i+ll/2)=snew(i)
  504. 210    continue
  505. c
  506. #ifdef SUNGRAPH
  507. c
  508. c *save high-passed future input in file 'ifile_ofile'
  509.    call save_sg(ifile_hp_vid, ssub, ll,'save ifile_hp_vid')
  510. #endif
  511. c
  512. c
  513. c    *** LPC spectral analysis (open loop)
  514. c    NOTE:  Autocorrelation was found superior to covariance analysis!
  515. c
  516.    if (anal .eq. 'autohf') then
  517.       call autohf(snew, hamw, ll, no, lambda, omega, fcn, rcn)
  518.    end if
  519. c
  520. #ifdef SUNGRAPH
  521. c *save rc's in file 'rc'
  522.   do 220 i=1,no
  523.      call save_sg(rc_vid(i), rcn(i), 1, 'save rc_vid')
  524. 220   continue
  525. #endif
  526. c
  527. c *pc -> lsp (new)
  528.    call pctolsp2(fcn, no, newfreq, lspflag)
  529.    if(lspflag) then
  530.       print *,' celp:  Bad "new" lsp at frame: ',frame
  531.       print 1030, 'lsp: ', (newfreq(i), i = 1, no)
  532.       print 1035, 'pc:  ', (fcn(i), i = 1, no+1)
  533.       print 1030, 'rc:  ', (rcn(i), i = 1, no)
  534. 1030       format(1x, a, 10f9.5)
  535. 1035       format(1x, a, 11f9.5)
  536.    end if
  537. c
  538. c *save unquantized lsp
  539.    do 230 i=1,no
  540.       unqfreq(i)=newfreq(i)
  541. 230    continue
  542. cc *quantize lsp's
  543.    if (stype .ne. 'none') then
  544.       if (stype .eq. 'kang') then
  545.  call lsp34(newfreq,no,sbits,findex)
  546.       end if
  547.    else
  548.       print *,'*** No lsp quantization'
  549.    end if
  550. c
  551. #ifdef SUNGRAPH
  552. c *save future lsp variables in file 'lsp1'
  553.    do 233 i=1,no
  554.       call save_sg(lsp_vid(i), unqfreq(i), 1,'save lsp_vid')
  555. 233    continue
  556. c
  557. c *save future qlsp variables in file 'lsp1'
  558.    do 237 i=1,no
  559.       call save_sg(qlsp_vid(i), newfreq(i), 1,'save qlsp_vid')
  560. 237    continue
  561. c
  562. #endif
  563. c *measure spectral distortion
  564. c *UNNECESSARY, used for diagnostics
  565.    call specdist(unqfreq, newfreq, dm2, sumdm2, framedm2)
  566. c
  567. c *pack lsp indices in bit stream array
  568.    do 240 i=1,no
  569.       call pack(findex(i),sbits(i),streambits,stream,pointer)
  570. 240    continue
  571. c
  572. c *linearly interpolate LSPs for each subframe
  573.    call intanaly(newfreq, nn, lsp)
  574. c
  575. c    *** for each subframe, search stochastic & adaptive code books
  576. c
  577.    k = 1
  578.    do 300 i = 1, nn
  579. #ifdef SUNGRAPH
  580. c *save interpolated lsp's in file 'lsp2'
  581.       do 310 j=1,no
  582.          call save_sg(lsp_analy_vid(j), lsp(j,i), 1,'save lsp_analy_vid')
  583. 310       continue
  584. #endif
  585.       call lsptopc(lsp(1,i),fci)
  586.       do 320 j=1,no+1
  587.          fc(j)=fci(j)
  588. 320       continue
  589.               nseg = nseg + 1
  590. c
  591. c       *** code book & pitch searches
  592. c
  593.               call csub(ssub(k),v(k),l,lp)
  595. #ifdef SUNGRAPH
  596. c *save code book index in file 'cbindex'
  597.       call save_sg(cb_index_vid, cbindex, 1,'save cb_index_vid')
  598. c
  599. c *save tau in file 'pitch'
  600.       call save_sg(pitch_tau_vid,bb(1),1,'save pitch_tau_vid')
  601. #endif
  602. c
  603. c *pitch quantization tau
  604. c
  605. c *pack parameter indices in bit stream array
  606.       if (mod(i,2) .ne. 0) then
  607.                  call packtau(tauptr-minptr,pbits(1),streambits,
  608.      &               pdencode,stream,pointer)
  609.       else
  610.                  call pack(tauptr-minptr,pbits(2),streambits,
  611.      &            stream,pointer)
  612.       end if
  613.               call pack(pindex,pbits(3),streambits,stream,pointer)
  614.       cbindex=cbindex-1
  615.       call pack(cbindex,cbbits,streambits,stream,pointer)
  616.       call pack(gindex,cbgbits,streambits,stream,pointer)
  617. c
  618. c *decode parameters for analysis by synthesis
  619.       cbindex=cbindex+1
  620. c
  621. c  *pitch synthesis (UNNECESSARY, used for diagnostics)
  622.       call pitchvq(v(k),l,dpa,idb,bb,'long')
  623. c
  624. c *pitch pre filter (UNNECESSARY, used for diagnostics)
  625.       if (prewt.ne.0.0) call prefilt(v(k),l,dppa)
  626. c
  627. c *lpc synthesis (UNNECESSARY, used for diagnostics)
  628.       call polefilt(fci,no,dsa,v(k),l)
  629. c
  630.               k = k + l
  631. 300    continue
  632. c
  633. c    *** bit error protection
  634. c                       *extract bits to protect from stream array
  635.            if (protect) then
  636.               do 410 i=1,codelength2
  637.                  codeword(i)=stream(bitprotect(i))
  638. 410           continue
  639. c
  640. c                       *Hamming encode
  641.              call encodeham(codelength1,codelength2,hmatrix,
  642.      &                    paritybit,codeword)
  643. c
  644. c *pack future bit
  645.       call pack(0,1,streambits,stream,pointer)
  646. c
  647. c                       *pack parity bits
  648.               do 420 i=1,paritylength
  649.                  call pack(codeword(codelength2+i),1,streambits,
  650.      &                     stream,pointer)
  651. 420           continue
  652. cnotused  call pack(paritybit,1,streambits,stream,pointer)
  653. c
  654. c *toggle and pack sync bit
  655.       sync = sync .xor. 1
  656.       call pack(sync,1,streambits,stream,pointer)
  657.            end if
  658. #ifdef SUNGRAPH
  659. c *save stream array in channel file
  660.    call save_sg(channel_vid, stream, streambits,'save channel_vid')
  661. c
  662. c *save stream array in stream_error file
  663.    call save_sg(stream_vid, stream, streambits,'save stream_vid')
  664. #endif
  665. c
  666. c *save stream
  667.    do 430 i=1,streambits
  668.       savestream(i)=stream(i)
  669. 430    continue
  670. c
  671. c *permute bitstream
  672.    do 440 i=1,streambits
  673.       pstream(i)=stream(bitpermute(i))
  674. 440    continue       
  675. c
  676. c *save stream in Dave's format
  677.    call puthex(streambits,pstream,line)
  678.    write (25,1040) line
  679. 1040    format (a36)
  680. c
  681. c
  682. c *** synthesizer only <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  683. c
  684. else
  685. c
  686. c    *** CHANNEL ..................................................
  687. c
  688. c     *read in channel file (if synthesis only)
  689.    read (25,1050,end=444) line
  690. 1050    format(a)
  691.    call gethex (streambits,pstream,line)
  692. c
  693. c *** end analyzer/synthesiser only differences >>>>>>>>>>>>>>>>>>>
  694. c
  695. end if
  696. c
  697. c *unpermute bitstream
  698. do 450 i=1,streambits
  699.    stream(bitpermute(i))=pstream(i)
  700. 450 continue
  701. c
  702. c    *corrupt the channel with bit errors 
  703. c
  704. call biterror(ber,mask,frame,stream,streambits,error,total)
  705. #ifdef SUNGRAPH
  706. c *save stream array in stream_error file
  707. if (analy) call save_sg(stream_error_vid, stream, streambits,
  708.      &                       'save stream_error_vid')
  709. #endif
  710. c
  711. c    *** SYNTHESIS ..................................................
  712. c
  713. c                       *unpack parity bits
  714. if (protect) then
  715.    pointer=pointer-paritylength-2
  716.    if (.not. analy) pointer=139
  717.    do 510 i=1,paritylength
  718.       call unpack(stream,1,streambits,codeword(codelength2+i),pointer)
  719. 510        continue
  720. c          call unpack(stream,1,streambits,paritybit,pointer)
  721. c
  722. c *extract code word from stream array
  723.            do 520 i=1,codelength2
  724.               codeword(i)=stream(bitprotect(i))
  725. 520        continue
  726. c
  727. c *repack future bit (remains constant for now)
  728.    codeword(11)=0
  729. c
  730. c                       *Hamming decode
  731.            call decodeham(codelength1,codelength2,hmatrix,
  732.      &        syndrometable,paritybit,codeword,twoerror,syndrome)
  733. c
  734. c *disable parity check (if parity not used)
  735.    twoerror = .false.
  736. c
  737. c *bit error rate estimator (running avg of bad syndromes)
  738.    if (syndrome .ne. 0) syndrome = 1
  739.    syndavg = (1.-(1./syndrun))*syndavg+(1./syndrun)*float(syndrome)
  740. c
  741. c                       *repack protected bits
  742.            do 530 i=1,codelength2
  743.               stream(bitprotect(i))=codeword(i)
  744. 530        continue
  745. #ifdef SUNGRAPH
  746. c *save stream array in stream_error file
  747.    if (analy) call save_sg(stream_error_s_vid, stream, streambits,
  748.      &                         'save stream_error_s_vid')
  749. #endif
  750. c
  751. c *compare saved stream with channel stream
  752.    errorcount=0
  753.    do 540 i=1,streambits
  754.       if (savestream(i) .ne. stream(i)) then
  755.  iint=1
  756.  errorcount=errorcount+1
  757. #ifdef SUNGRAPH
  758.  if (analy) call save_sg(bit_error_vid,iint,1,
  759.      &                  'save bit_error_vid')
  760. #endif
  761.       else
  762.  iint=0
  763. #ifdef SUNGRAPH
  764.     if (analy) call save_sg(bit_error_vid,iint,1,
  765.      &                       'save bit_error_vid')
  766. #endif
  767.       end if
  768. 540    continue
  769.    perrorcount=0
  770.    do 550 i=1,codelength2
  771.       if (savestream(bitprotect(i)) .ne. stream(bitprotect(i)))
  772.      &               perrorcount=perrorcount+1
  773. 550    continue
  774. c
  775. c *frame repeat if two errors detected in code word
  776. copt           if (twoerror) then
  777. copt       print *,' celp:  two errors have occured in frame',frame
  778. copt             do 560 i=1,streambits
  779. copt                stream(i)=oldstream(i)
  780. copt560          continue
  781. copt             twoerror = .false.
  782. copt           end if
  783.            do 570 i=1,streambits
  784.               oldstream(i)=stream(i)
  785. 570        continue
  786. end if
  787. c
  788. pointer=0
  789. c
  790. c *unpack data stream
  791. do 580 i=1,no
  792.    call unpack(stream,sbits(i),streambits,findex(i),pointer)
  793. 580 continue
  794. c *decode lsp's
  795. if (stype .ne. 'none') then
  796.    if (stype .eq. 'kang') then
  797.       call lspdecode34(findex,no,newfreq)
  798.    else
  799.       stop 'Problem with lspdecoding'
  800.    end if
  801. end if
  802. c *interpolate spectrum lsp's for nn subframes
  803. call intsynth(newfreq, nn, lsp, twoerror, syndavg, predgain,
  804.      &    preddiff, predavg, mindiff)
  805. #ifdef SUNGRAPH
  806. c *save spectrum smoothing parameters in file 'rc'
  807. if (analy) then
  808.    call save_sg(predgain_vid, predgain, 1, 'save predgain_vid')
  809.    call save_sg(preddiff_vid, preddiff, 1, 'save preddiff_vid')
  810.    call save_sg(predavg_vid, predavg, 1, 'save predavg_vid')
  811.    call save_sg(mindiff_vid, mindiff, 1, 'save mindiff_vid')
  812. end if
  813. #endif
  814. c
  815. c *decode all code book and pitch parameters
  816. bitpointer=pointer
  817. call dcodtau(pbits(1),pbits(2),bitsum1,bitsum2,bitpointer,nn,
  818.      &    stream,streambits,pddecode,pdtabi,taus)
  819. call dcodpg(pbits(3),ptype,bitsum1,bitsum2,bitpointer,nn,stream,
  820.      &    streambits,pgs)
  821. call dcodcbi(cbbits,bitsum1,bitsum2,bitpointer,nn,stream,
  822.      &    streambits,cbi)
  823. call dcodcbg(cbgbits,cbgtype,bitsum1,bitsum2,bitpointer,nn,stream,
  824.      &    streambits,cbg)
  825. c
  826. c    *** synthesize each subframe
  827. c
  828. if (analy) nseg = nseg - nn
  829. k = 1
  830. do 600 i = 1, nn
  831.    nseg = nseg + 1
  832. c
  833. #ifdef SUNGRAPH
  834. c *save interpolated lsp's in file 'lsp2'
  835.    if (analy) then
  836.    do 610 j=1,no
  837.       call save_sg(lsp_synth_vid(j), lsp(j,i), 1,'save lsp_synth_vid')
  838. 610    continue
  839.    end if
  840. #endif
  841. c
  842. c *decode values for subframe 
  843.    cbindex = cbi(i)
  844.    decodedgain = cbg(i)
  845.    if (protect) call smoothcbgain(decodedgain,twoerror,syndavg,
  846.      & cbg,i)
  847. c
  848. c *code book synthesis
  849.    call vdecode(decodedgain,l,vdecoded(k))
  850. c
  851.            if (protect) call smoothtau(taus(i),twoerror,
  852.      &                                 syndavg,taus(3),i)
  853.    bb(1)=taus(i)
  854.    bb(3)=pgs(i)
  855.    if (protect) call smoothpgain(bb(3),twoerror,syndavg,pgs,i)
  856. c
  857. #ifdef SUNGRAPH
  858. c *save synthesis parameters in sungraph files
  859.    if (analy) call save_sg(cb_index_synth_vid, cbindex, 1,
  860.      &       'save cb_index_synth_vid')
  861.    if (analy) call save_sg(pitch_tau_synth_vid, bb(1), 1,
  862.      &       'save pitch_tau_synth_vid')
  863.    if (analy) call save_sg(cb_qgain_synth_vid, decodedgain, 1,
  864.      &       'save cb_qgain_synth_vid')
  865.    if (analy) call save_sg(pitch_qgain_synth_vid, bb(3), 1, 
  866.      &       'save pitch_qgain_synth_vid')
  867. #endif
  868. c
  869. c *pitch synthesis
  870.    call pitchvq(vdecoded(k),l,dps,idb,bb,'long')
  871. c
  872. c *pitch pre filter (synthesis)
  873.    if (prewt.ne.0.0) call prefilt(vdecoded(k),l,dpps)
  874. c
  875. c *convert lsp's to pc's
  876.    call lsptopc(lsp(1,i),fci)
  877. c
  878. c *lpc synthesis
  879.    call polefilt(fci,no,dss,vdecoded(k),l)
  880. c
  881. c       *** check analysis versus synthesis speech
  882. c
  883.    flag = .false.
  884.    do 620 j = 0, l-1
  885.       if (abs(v(k+j)-vdecoded(k+j)) .gt. 1.e-6) flag=.true.
  886. 620    continue
  887.    if (flag .and. analy) print *,
  888.      &               ' celp:  Speech mismatch at frame',frame
  889. c
  890. c       *** write nonpostfiltered output speech disk files
  891. c
  892.    do 630 j = 0, l-1 
  893.       vdecoded(k+j) = descale * vdecoded(k+j)
  894.       npf(k+j)=nint(max(-32768.,min (32767., vdecoded(k+j))))
  895. canaly       npf(k+j)=nint(max(-32768.,min (32767., v(k+j))))
  896. 630    continue
  897. c *write npf output speech file "ofilenpf"
  898. cnpf    if (iodisk(2,2,ofile(1:oflen)//'npf.spd',nrec2,npf(k),l).ne.l)
  899. cnpf     &          go to 888
  900. #ifdef SUNGRAPH
  901. c *sungraph's npf output speech file 'ifile_ofile'
  902.    if (analy) call save_sg(ofile_npf_vid,npf(k),l,
  903.      & 'save ofile_npf_vid')
  904. #endif
  905. c
  906. c       *** calculate the average segmental SNR
  907. c *UNNECESSARY, used for diagnostics
  908. c
  909.    if (analy) call segsnr(ssub(k),npf(k),l,sumsnr,framesnr,snrflag)
  910. c
  911. c       *** calculate distortions/distances (log spectral error, etc.)
  912. c *UNNECESSARY, used for diagnostics
  913. c
  914.    if (snrflag .and. analy) call distortion(ssub(k),npf(k),hamws,
  915.      &                    l,no,dm,sumdm,framedm)
  916. c
  917. c       *** post filtering
  918. c
  919.   call postfilt(vdecoded(k),l,alpha,beta,ip,op,dp1,dp2,dp3,dp4)
  920. canaly   call postfilt(v(k),l,alpha,beta,ip,op,dp1,dp2,dp3,dp4)
  921. clocomp   call postfilt2(vdecoded(k),l,alpha,beta,ip,op,dp1,dp2,dp3,dp4)
  922. c
  923. c       *** test for output speech clipping
  924. c
  925. 636    if (clip(vdecoded(k), l)) then
  926. c *frame repeat & recall synthesizer?
  927. c *or scale vdecoded?
  928.       print *,' celp:  Clipping detected at frame', frame
  929.       do 635 j = 0, l-1
  930.  vdecoded(k+j) = 0.05*vdecoded(k+j)
  931. 635       continue
  932.       go to 636
  933.    end if
  934. c
  935. c       *** write postfiltered output speech disk files
  936. c
  937.   do 640 j = 0, l-1
  938.       pf(k+j)=nint(max(-32768.,min (32767., vdecoded(k+j))))
  939. canaly       pf(k+j)=nint(max(-32768.,min (32767., v(k+j))))
  940. 640    continue
  941. c
  942. c *write output speech file "ofile"
  943.    if (iodisk(2,1,ofile(1:oflen)//'.spd',nrec1,pf(k),l).ne.l)
  944.      &          go to 888
  945. #ifdef SUNGRAPH
  946. c *sungraph's output speech file 'ifile_ofile'
  947.    if (analy) call save_sg(ofile_pf_vid,pf(k),l,'save ofile_pf_vid')
  948. #endif
  949. c
  950. c       *** high pass filter output speech
  951. c
  952.    call zerofilt(ahpfo,2,dhpf1o,vdecoded(k),l)
  953.    call polefilt(bhpfo,2,dhpf2o,vdecoded(k),l)
  954. c
  955.    do 641 j = 0, l-1
  956.       pf(k+j)=nint(max(-32768.,min (32767., vdecoded(k+j))))
  957. 641    continue
  958. c
  959. c *write output speech file "ofilehpf"
  960.    if (iodisk(2,3,ofile(1:oflen)//'hpf.spd',nrec3,pf(k),l).ne.l)
  961.      &          go to 888
  962. c
  963. #ifdef SUNGRAPH
  964. c *sungraph's output speech file 'ifile_ofile'
  965.    if (analy) call save_sg(ofile_hpf_vid,pf(k),l,'save ofile_hpf_vid')
  966. #endif
  967. c
  968.    k = k + l
  969. 600 continue
  970. c    . . . . . . . . . . . . . . . . . . . . . . .
  971. c
  972. c *** end block
  973. c
  974. frame = frame + 1
  975. c
  976. c *display a propeller (rotating bar) once per frame
  977. call mark(0)
  978. cframe print *,'frame= ',frame
  979. if (.not. analy) go to 69
  980. c
  981. #ifdef SUNGRAPH
  982. include 'sgeblk.com'
  983. convex #include "sgeblk.com"
  984. #endif
  985. c
  986. c    *** shift new speech buffer into old speech buffer
  987. c
  988. c sold    snew
  989. c |-------------------|-------------------| snew
  990. c   |-------------------| 
  991. c    ssub
  992. c
  993. do 730 i = 1, ll
  994.    sold(i) = snew(i)
  995. 730 continue
  996. c
  997. c    *** frame finished, end loop
  998. go to 69
  999. c
  1000. c................ e n d  m a i n  l o o p ......(fortran 77 do while)......
  1001. c
  1002. c *** WRAP-UP ..................................................
  1003. c
  1004. c
  1005. c *** exit through here if a file write error occurs
  1006. c
  1007. 888 print *,' celp:  Error writing output file'
  1008. c
  1009. c *** finished reading & processing "ifile"
  1010. c
  1011. 999     continue
  1012. c
  1013. c *** closing cue and close files
  1014. c
  1015. do 810 i = 1, streambits
  1016.    stream(i) = 23456
  1017. 810 continue
  1018.         ierr=iodisk(-1,1,ofile(1:oflen)//'.spd',nrec1,pf,l)
  1019. cnpf ierr=iodisk(-1,2,ofile(1:oflen)//'npf.spd',nrec2,npf,l)
  1020. ierr=iodisk(-1,3,ofile(1:oflen)//'hpf.spd',nrec3,pf,l)
  1021. cc close (20)
  1022. cc close (21)
  1023. close (25)
  1024. c
  1025. #ifdef SUNGRAPH
  1026. include 'sgclose.com'
  1027. convex #include "sgclose.com"
  1028. #else
  1029. ierr=iodisk(-1,9,ifile,nrec9,iarf,ll)
  1030. #endif
  1031. if (total .ne. 0) then
  1032.    realerror=100.*(float(error)/float(total))
  1033. else
  1034.    realerror=0.0
  1035. end if
  1036. 444 call cliend(sumsnr,framesnr,realerror,
  1037.      &     sumdm,framedm,sumdm2,framedm2, analy)
  1038. if (.not. analy) print *,'Synthesis run successful'
  1039. stop
  1040. end