1#!/usr/bin/python
2
3# Audio Tools, a module and set of tools for manipulating audio data
4# Copyright (C) 2007-2014  Brian Langenberger
5
6# This program is free software; you can redistribute it and/or modify
7# it under the terms of the GNU General Public License as published by
8# the Free Software Foundation; either version 2 of the License, or
9# (at your option) any later version.
10
11# This program is distributed in the hope that it will be useful,
12# but WITHOUT ANY WARRANTY; without even the implied warranty of
13# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14# GNU General Public License for more details.
15
16# You should have received a copy of the GNU General Public License
17# along with this program; if not, write to the Free Software
18# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
19
20from audiotools.bitstream import BitstreamWriter
21from audiotools.bitstream import BitstreamRecorder
22from audiotools import BufferedPCMReader
23from hashlib import md5
24
25
26class STREAMINFO(object):
27    def __init__(self, minimum_block_size, maximum_block_size,
28                 minimum_frame_size, maximum_frame_size,
29                 sample_rate, channels, bits_per_sample,
30                 total_pcm_frames, md5sum):
31        self.minimum_block_size = minimum_block_size
32        self.maximum_block_size = maximum_block_size
33        self.minimum_frame_size = minimum_frame_size
34        self.maximum_frame_size = maximum_frame_size
35        self.sample_rate = sample_rate
36        self.channels = channels
37        self.bits_per_sample = bits_per_sample
38        self.total_pcm_frames = total_pcm_frames
39        self.md5sum = md5sum
40
41    def write(self, writer):
42        writer.build("16u 16u 24u 24u 20u 3u 5u 36U 16b",
43                     [self.minimum_block_size,
44                      self.maximum_block_size,
45                      self.minimum_frame_size,
46                      self.maximum_frame_size,
47                      self.sample_rate,
48                      self.channels - 1,
49                      self.bits_per_sample - 1,
50                      self.total_pcm_frames,
51                      self.md5sum.digest()])
52
53    def input_update(self, framelist):
54        self.total_pcm_frames += framelist.frames
55        self.md5sum.update(framelist.to_bytes(False, True))
56
57    def output_update(self, flac_frame):
58        self.minimum_frame_size = min(self.minimum_frame_size,
59                                      flac_frame.bytes())
60        self.maximum_frame_size = max(self.maximum_frame_size,
61                                      flac_frame.bytes())
62
63
64class Encoding_Options(object):
65    def __init__(self, block_size=4096, max_lpc_order=8,
66                 adaptive_mid_side=False, mid_side=True,
67                 exhaustive_model_search=False,
68                 max_residual_partition_order=5,
69                 max_rice_parameter=14):
70        self.block_size = block_size
71        self.max_lpc_order = max_lpc_order
72        self.adaptive_mid_side = adaptive_mid_side
73        self.mid_side = mid_side
74        self.exhaustive_model_search = exhaustive_model_search
75        self.max_residual_partition_order = max_residual_partition_order
76        self.max_rice_parameter = max_rice_parameter
77
78        if block_size <= 192:
79            self.qlp_precision = 7
80        elif block_size <= 384:
81            self.qlp_precision = 8
82        elif block_size <= 576:
83            self.qlp_precision = 9
84        elif block_size <= 1152:
85            self.qlp_precision = 10
86        elif block_size <= 2304:
87            self.qlp_precision = 11
88        elif block_size <= 4608:
89            self.qlp_precision = 12
90        else:
91            self.qlp_precision = 13
92
93
94def encode_flac(filename,
95                pcmreader,
96                block_size=4096,
97                max_lpc_order=8,
98                min_residual_partition_order=0,
99                max_residual_partition_order=5,
100                mid_side=True,
101                adaptive_mid_side=False,
102                exhaustive_model_search=False,
103                disable_verbatim_subframes=False,
104                disable_constant_subframes=False,
105                disable_fixed_subframes=False,
106                disable_lpc_subframes=False,
107                padding_size=4096):
108
109    current_offset = 0
110    frame_offsets = []
111
112    options = Encoding_Options(block_size,
113                               max_lpc_order,
114                               adaptive_mid_side,
115                               mid_side,
116                               exhaustive_model_search,
117                               max_residual_partition_order,
118                               14 if pcmreader.bits_per_sample <= 16 else 30)
119
120    streaminfo = STREAMINFO(block_size,
121                            block_size,
122                            (2 ** 24) - 1,
123                            0,
124                            pcmreader.sample_rate,
125                            pcmreader.channels,
126                            pcmreader.bits_per_sample,
127                            0, md5())
128
129    pcmreader = BufferedPCMReader(pcmreader)
130    output_file = open(filename, "wb")
131    writer = BitstreamWriter(output_file, False)
132
133    # write placeholder metadata blocks such as STREAMINFO and PADDING
134    writer.write_bytes("fLaC")
135    writer.build("1u 7u 24u", [0, 0, 34])
136    streaminfo_start = writer.getpos()
137    streaminfo.write(writer)
138
139    writer.build("1u 7u 24u", [1, 1, padding_size])
140    writer.write_bytes(b"\x00" * padding_size)
141
142    # walk through PCM reader's FrameLists
143    frame_number = 0
144    frame = pcmreader.read(block_size)
145
146    flac_frame = BitstreamRecorder(0)
147
148    while len(frame) > 0:
149        frame_offsets.append((current_offset, frame.frames))
150        streaminfo.input_update(frame)
151
152        flac_frame.reset()
153        encode_flac_frame(flac_frame, pcmreader, options, frame_number, frame)
154        current_offset += flac_frame.bytes()
155        streaminfo.output_update(flac_frame)
156
157        flac_frame.copy(writer)
158
159        frame_number += 1
160        frame = pcmreader.read(block_size)
161
162    # return to beginning of file and rewrite STREAMINFO block
163    writer.setpos(streaminfo_start)
164    streaminfo.write(writer)
165    writer.flush()
166    writer.close()
167
168    return frame_offsets
169
170
171def encode_flac_frame(writer, pcmreader, options, frame_number, frame):
172    crc16 = CRC16()
173    writer.add_callback(crc16.update)
174
175    if ((pcmreader.channels == 2) and (options.adaptive_mid_side or
176                                       options.mid_side)):
177        # calculate average/difference
178        average = [(c0 + c1) // 2 for c0, c1 in zip(frame.channel(0),
179                                                    frame.channel(1))]
180        difference = [c0 - c1 for c0, c1 in zip(frame.channel(0),
181                                                frame.channel(1))]
182
183        # try different subframes based on encoding options
184        left_subframe = BitstreamRecorder(0)
185        encode_subframe(left_subframe, options,
186                        pcmreader.bits_per_sample, list(frame.channel(0)))
187
188        right_subframe = BitstreamRecorder(0)
189        encode_subframe(right_subframe, options,
190                        pcmreader.bits_per_sample, list(frame.channel(1)))
191
192        average_subframe = BitstreamRecorder(0)
193        encode_subframe(average_subframe, options,
194                        pcmreader.bits_per_sample, average)
195
196        difference_subframe = BitstreamRecorder(0)
197        encode_subframe(difference_subframe, options,
198                        pcmreader.bits_per_sample + 1, difference)
199
200        # write best header/subframes to disk
201        if options.mid_side:
202            if ((left_subframe.bits() + right_subframe.bits()) <
203                min(left_subframe.bits() + difference_subframe.bits(),
204                    difference_subframe.bits() + right_subframe.bits(),
205                    average_subframe.bits() + difference_subframe.bits())):
206                write_frame_header(writer, pcmreader, frame_number, frame, 0x1)
207                left_subframe.copy(writer)
208                right_subframe.copy(writer)
209            elif (left_subframe.bits() <
210                  min(right_subframe.bits(), difference_subframe.bits())):
211                write_frame_header(writer, pcmreader, frame_number, frame, 0x8)
212                left_subframe.copy(writer)
213                difference_subframe.copy(writer)
214            elif right_subframe.bits() < average_subframe.bits():
215                write_frame_header(writer, pcmreader, frame_number, frame, 0x9)
216                difference_subframe.copy(writer)
217                right_subframe.copy(writer)
218            else:
219                write_frame_header(writer, pcmreader, frame_number, frame, 0xA)
220                average_subframe.copy(writer)
221                difference_subframe.copy(writer)
222        else:
223            if (((left_subframe.bits() + right_subframe.bits()) <
224                 (average_subframe.bits() + difference_subframe.bits()))):
225                write_frame_header(writer, pcmreader, frame_number, frame, 0x1)
226                left_subframe.copy(writer)
227                right_subframe.copy(writer)
228            else:
229                write_frame_header(writer, pcmreader, frame_number, frame, 0xA)
230                average_subframe.copy(writer)
231                difference_subframe.copy(writer)
232    else:
233        write_frame_header(writer, pcmreader, frame_number, frame,
234                           pcmreader.channels - 1)
235
236        for i in range(frame.channels):
237            encode_subframe(writer,
238                            options,
239                            pcmreader.bits_per_sample,
240                            list(frame.channel(i)))
241
242    writer.byte_align()
243    writer.pop_callback()
244    writer.write(16, int(crc16))
245
246
247def write_frame_header(writer, pcmreader, frame_number, frame,
248                       channel_assignment):
249    crc8 = CRC8()
250    writer.add_callback(crc8.update)
251    writer.write(14, 0x3FFE)
252    writer.write(1, 0)
253    writer.write(1, 0)
254    encoded_block_size = {192: 1,
255                          256: 8,
256                          512: 9,
257                          576: 2,
258                          1024: 10,
259                          1152: 3,
260                          2048: 11,
261                          2304: 4,
262                          4096: 12,
263                          4608: 5,
264                          8192: 13,
265                          16384: 14,
266                          32768: 15}.get(frame.frames, None)
267    if encoded_block_size is None:
268        if frame.frames <= 256:
269            encoded_block_size = 6
270        elif frame.frames <= 65536:
271            encoded_block_size = 7
272        else:
273            encoded_block_size = 0
274    writer.write(4, encoded_block_size)
275
276    encoded_sample_rate = {8000: 4,
277                           16000: 5,
278                           22050: 6,
279                           24000: 7,
280                           32000: 8,
281                           44100: 9,
282                           48000: 10,
283                           88200: 1,
284                           96000: 11,
285                           176400: 2,
286                           192000: 3}.get(pcmreader.sample_rate, None)
287    if encoded_sample_rate is None:
288        if ((((pcmreader.sample_rate % 1000) == 0) and
289             (pcmreader.sample_rate <= 255000))):
290            encoded_sample_rate = 12
291        elif (((pcmreader.sample_rate % 10) == 0) and
292              (pcmreader.sample_rate <= 655350)):
293            encoded_sample_rate = 14
294        elif pcmreader.sample_rate <= 65535:
295            encoded_sample_rate = 13
296        else:
297            encoded_sample_rate = 0
298    writer.write(4, encoded_sample_rate)
299
300    writer.write(4, channel_assignment)
301
302    encoded_bps = {8: 1,
303                   12: 2,
304                   16: 4,
305                   20: 5,
306                   24: 6}.get(pcmreader.bits_per_sample, 0)
307    writer.write(3, encoded_bps)
308
309    writer.write(1, 0)
310    write_utf8(writer, frame_number)
311
312    if encoded_block_size == 6:
313        writer.write(8, frame.frames - 1)
314    elif encoded_block_size == 7:
315        writer.write(16, frame.frames - 1)
316
317    if encoded_sample_rate == 12:
318        writer.write(8, pcmreader.sample_rate % 1000)
319    elif encoded_sample_rate == 13:
320        writer.write(16, pcmreader.sample_rate)
321    elif encoded_sample_rate == 14:
322        writer.write(16, pcmreader.sample_rate % 10)
323
324    writer.pop_callback()
325    writer.write(8, int(crc8))
326
327
328def write_utf8(writer, value):
329    if value <= 127:
330        writer.write(8, value)
331    else:
332        if value <= 2047:
333            total_bytes = 2
334        elif value <= 65535:
335            total_bytes = 3
336        elif value <= 2097151:
337            total_bytes = 4
338        elif value <= 67108863:
339            total_bytes = 5
340        elif value <= 2147483647:
341            total_bytes = 6
342        else:
343            raise ValueError("UTF-8 value too large")
344
345        shift = (total_bytes - 1) * 6
346        writer.unary(0, total_bytes)
347        writer.write(7 - total_bytes, value >> shift)
348        shift -= 6
349        while shift >= 0:
350            writer.write(2, 2)
351            writer.write(6, (value >> shift) & 0x3F)
352            shift -= 6
353
354
355def encode_subframe(writer, options, bits_per_sample, samples):
356    def all_identical(l):
357        if len(l) == 1:
358            return True
359        else:
360            for i in l[1:]:
361                if i != l[0]:
362                    return False
363            else:
364                return True
365
366    def wasted(s):
367        w = 0
368        while (s & 1) == 0:
369            w += 1
370            s >>= 1
371        return w
372
373    if all_identical(samples):
374        encode_constant_subframe(writer, bits_per_sample, samples[0])
375    else:
376        # account for wasted BPS, if any
377        wasted_bps = 2 ** 32
378        for sample in samples:
379            if sample != 0:
380                wasted_bps = min(wasted_bps, wasted(sample))
381                if wasted_bps == 0:
382                    break
383
384        if wasted_bps == 2 ** 32:
385            # all samples are 0
386            wasted_bps = 0
387        elif wasted_bps > 0:
388            samples = [s >> wasted_bps for s in samples]
389
390        fixed_subframe = BitstreamRecorder(0)
391        encode_fixed_subframe(fixed_subframe,
392                              options,
393                              wasted_bps,
394                              bits_per_sample,
395                              samples)
396
397        if options.max_lpc_order > 0:
398            lpc_subframe = BitstreamRecorder(0)
399            encode_lpc_subframe(lpc_subframe,
400                                options,
401                                wasted_bps,
402                                bits_per_sample,
403                                samples)
404
405            if (((bits_per_sample * len(samples)) <
406                 min(fixed_subframe.bits(), lpc_subframe.bits()))):
407                encode_verbatim_subframe(writer, wasted_bps,
408                                         bits_per_sample, samples)
409            elif fixed_subframe.bits() < lpc_subframe.bits():
410                fixed_subframe.copy(writer)
411            else:
412                lpc_subframe.copy(writer)
413        else:
414            if (bits_per_sample * len(samples)) < fixed_subframe.bits():
415                encode_verbatim_subframe(writer, wasted_bps,
416                                         bits_per_sample, samples)
417            else:
418                fixed_subframe.copy(writer)
419
420
421def encode_constant_subframe(writer, bits_per_sample, sample):
422    # write frame header
423    writer.build("1p 6u 1u", [0, 0])
424
425    # write frame data
426    writer.write_signed(bits_per_sample, sample)
427
428
429def encode_verbatim_subframe(writer, wasted_bps, bits_per_sample, samples):
430    # write frame header
431    writer.build("1p 6u", [1])
432    if wasted_bps > 0:
433        writer.write(1, 1)
434        writer.unary(1, wasted_bps - 1)
435    else:
436        writer.write(1, 0)
437
438    # write frame data
439    writer.build(("%ds" % (bits_per_sample - wasted_bps)) * len(samples),
440                 samples)
441
442
443def encode_fixed_subframe(writer, options, wasted_bps, bits_per_sample,
444                          samples):
445    def next_order(residuals):
446        return [(x - y) for (x, y) in zip(residuals[1:], residuals)]
447
448    # decide which subframe order to use
449    residuals = [samples]
450    total_error = [sum(map(abs, residuals[-1][4:]))]
451
452    if len(samples) > 4:
453        for order in range(1, 5):
454            residuals.append(next_order(residuals[-1]))
455            total_error.append(sum(map(abs, residuals[-1][4 - order:])))
456
457        for order in range(4):
458            if total_error[order] < min(total_error[order + 1:]):
459                break
460        else:
461            order = 4
462    else:
463        order = 0
464
465    # then write the subframe to disk
466
467    # write subframe header
468    writer.build("1p 3u 3u", [1, order])
469    if wasted_bps > 0:
470        writer.write(1, 1)
471        writer.unary(1, wasted_bps - 1)
472    else:
473        writer.write(1, 0)
474
475    # write warm-up samples
476    for sample in samples[0:order]:
477        writer.write_signed(bits_per_sample - wasted_bps, sample)
478
479    # write residual block
480    encode_residuals(writer, options, order, len(samples), residuals[order])
481
482
483def encode_residuals(writer, options, order, block_size, residuals):
484    # first, determine the best set of residual partitions to use
485    best_porder = 0
486    best_size = 2 ** 31
487
488    for porder in range(0, options.max_residual_partition_order + 1):
489        if (block_size % (2 ** porder)) == 0:
490            unencoded_residuals = residuals[:]
491            partitions = []
492            for p in range(0, 2 ** porder):
493                if p == 0:
494                    partition_size = block_size // (2 ** porder) - order
495                else:
496                    partition_size = block_size // (2 ** porder)
497                partitions.append(unencoded_residuals[0:partition_size])
498                unencoded_residuals = unencoded_residuals[partition_size:]
499
500            rice_parameters = [best_rice_parameter(options, p)
501                               for p in partitions]
502
503            encoded_partitions = [encode_residual_partition(r, p)
504                                  for (r, p) in zip(rice_parameters,
505                                                    partitions)]
506
507            partition_bit_size = sum([4 + p.bits()
508                                      for p in encoded_partitions])
509
510            if partition_bit_size < best_size:
511                best_porder = porder
512                best_size = partition_bit_size
513                best_parameters = rice_parameters
514                best_encoded_partitions = encoded_partitions
515
516    # then output those residual partitions into a single block
517    if max(best_parameters) > 14:
518        coding_method = 1
519    else:
520        coding_method = 0
521
522    writer.write(2, coding_method)
523    writer.write(4, best_porder)
524    for (rice, partition) in zip(best_parameters, best_encoded_partitions):
525        if coding_method == 0:
526            writer.write(4, rice)
527        else:
528            writer.write(5, rice)
529        partition.copy(writer)
530
531
532def best_rice_parameter(options, residuals):
533    partition_sum = sum(map(abs, residuals))
534    rice_parameter = 0
535    while (len(residuals) * (2 ** rice_parameter)) < partition_sum:
536        if rice_parameter < options.max_rice_parameter:
537            rice_parameter += 1
538        else:
539            return options.max_rice_parameter
540
541    return rice_parameter
542
543
544def encode_residual_partition(rice_parameter, residuals):
545    partition = BitstreamRecorder(0)
546    for residual in residuals:
547        if residual >= 0:
548            unsigned = residual << 1
549        else:
550            unsigned = ((-residual - 1) << 1) | 1
551        MSB = unsigned >> rice_parameter
552        LSB = unsigned - (MSB << rice_parameter)
553        partition.unary(1, MSB)
554        partition.write(rice_parameter, LSB)
555
556    return partition
557
558
559def encode_lpc_subframe(writer, options, wasted_bps, bits_per_sample,
560                        samples):
561    """computes the best LPC subframe from the given samples
562    according to the options and writes it to the given BitstreamWriter"""
563
564    # window signal
565    windowed = [(sample * tukey) for (sample, tukey) in
566                zip(samples, tukey_window(len(samples), 0.5))]
567
568    # compute autocorrelation values
569    if len(samples) > (options.max_lpc_order + 1):
570        autocorrelation_values = [
571            sum(x * y for x, y in zip(windowed, windowed[lag:]))
572            for lag in range(0, options.max_lpc_order + 1)]
573    else:
574        # not enough samples, so build LPC with dummy coeffs
575        write_lpc_subframe(writer=writer,
576                           options=options,
577                           wasted_bps=wasted_bps,
578                           bits_per_sample=bits_per_sample,
579                           order=1,
580                           qlp_precision=options.qlp_precision,
581                           qlp_shift_needed=0,
582                           qlp_coefficients=[0],
583                           samples=samples)
584        return
585
586    # compute LP coefficients from autocorrelation values
587    if ((len(autocorrelation_values) > 1) and
588        (set(autocorrelation_values) != {0.0})):
589        (lp_coefficients, error) = \
590            compute_lp_coefficients(autocorrelation_values)
591    else:
592        # all autocorrelation values are zero
593        # so build LPC with dummy coeffs
594        write_lpc_subframe(writer=writer,
595                           options=options,
596                           wasted_bps=wasted_bps,
597                           bits_per_sample=bits_per_sample,
598                           order=1,
599                           qlp_precision=options.qlp_precision,
600                           qlp_shift_needed=0,
601                           qlp_coefficients=[0],
602                           samples=samples)
603        return
604
605    if not options.exhaustive_model_search:
606        # if not performaing an exhaustive model search
607        # estimate which set of LP coefficients is best
608        # and use those to build subframe
609
610        order = estimate_best_lpc_order(options,
611                                        len(samples),
612                                        bits_per_sample,
613                                        error)
614        (qlp_coefficients, qlp_shift_needed) = \
615            quantize_coefficients(options.qlp_precision,
616                                  lp_coefficients,
617                                  order)
618
619        write_lpc_subframe(writer=writer,
620                           options=options,
621                           wasted_bps=wasted_bps,
622                           bits_per_sample=bits_per_sample,
623                           order=order,
624                           qlp_precision=options.qlp_precision,
625                           qlp_shift_needed=qlp_shift_needed,
626                           qlp_coefficients=qlp_coefficients,
627                           samples=samples)
628    else:
629        # otherwise, build all possible subframes
630        # and return the one which is actually the smallest
631
632        best_subframe_size = 2 ** 32
633        best_subframe = BitstreamRecorder(0)
634
635        for order in range(1, options.max_lpc_order + 1):
636            (qlp_coefficients, qlp_shift_needed) = \
637                quantize_coefficients(options.qlp_precision,
638                                      lp_coefficients,
639                                      order)
640
641            subframe = BitstreamRecorder(0)
642            write_lpc_subframe(writer=subframe,
643                               options=options,
644                               wasted_bps=wasted_bps,
645                               bits_per_sample=bits_per_sample,
646                               order=order,
647                               qlp_precision=options.qlp_precision,
648                               qlp_shift_needed=qlp_shift_needed,
649                               qlp_coefficients=qlp_coefficients,
650                               samples=samples)
651            if subframe.bits() < best_subframe_size:
652                best_subframe = subframe
653        else:
654            best_subframe.copy(writer)
655
656
657def tukey_window(sample_count, alpha):
658    from math import cos, pi
659
660    window1 = (alpha * (sample_count - 1)) // 2
661    window2 = (sample_count - 1) * (1 - (alpha // 2))
662
663    for n in range(0, sample_count):
664        if n <= window1:
665            yield (0.5 *
666                   (1 +
667                    cos(pi * (((2 * n) / (alpha * (sample_count - 1))) - 1))))
668        elif n <= window2:
669            yield 1.0
670        else:
671            yield (0.5 *
672                   (1 +
673                    cos(pi * (((2 * n) / (alpha * (sample_count - 1))) -
674                              (2 / alpha) + 1))))
675
676
677def compute_lp_coefficients(autocorrelation):
678    maximum_lpc_order = len(autocorrelation) - 1
679
680    k0 = autocorrelation[1] / autocorrelation[0]
681    lp_coefficients = [[k0]]
682    error = [autocorrelation[0] * (1 - k0 ** 2)]
683
684    for i in range(1, maximum_lpc_order):
685        ki = (autocorrelation[i + 1] -
686              sum([x * y for (x, y) in
687                   zip(lp_coefficients[i - 1],
688                       reversed(autocorrelation[1:i + 1]))])) / error[i - 1]
689
690        lp_coefficients.append([c1 - (ki * c2) for (c1, c2) in
691                                zip(lp_coefficients[i - 1],
692                                    reversed(lp_coefficients[i - 1]))] + [ki])
693        error.append(error[i - 1] * (1 - ki ** 2))
694
695    return (lp_coefficients, error)
696
697
698def estimate_best_lpc_order(options, block_size, bits_per_sample, error):
699    """returns an order integer of the best LPC order to use"""
700
701    from math import log
702
703    error_scale = log(2) ** 2
704    best_order = 0
705    best_subframe_bits = 1e32
706    for i in range(options.max_lpc_order):
707        order = i + 1
708        if error[i] > 0.0:
709            header_bits = order * (bits_per_sample + options.qlp_precision)
710            bits_per_residual = max((log(error[i] * error_scale) /
711                                     (log(2) * 2)), 0.0)
712            estimated_subframe_bits = (header_bits +
713                                       bits_per_residual *
714                                       (block_size - order))
715            if estimated_subframe_bits < best_subframe_bits:
716                best_order = order
717                best_subframe_bits = estimated_subframe_bits
718        elif error[i] == 0.0:
719            return order
720    else:
721        return best_order
722
723
724def quantize_coefficients(qlp_precision, lp_coefficients, order):
725    """returns a (qlp_coeffs, qlp_shift_needed) pair
726    where qlp_coeffs is a list of ints
727    and qlp_shift_needed is a non-negative integer"""
728
729    from math import log
730
731    l = max(map(abs, lp_coefficients[order - 1]))
732    if l > 0:
733        qlp_shift_needed = min((qlp_precision - 1) -
734                               (int(log(l) / log(2)) - 1) - 1,
735                               (2 ** 4) - 1)
736    else:
737        qlp_shift_needed = 0
738    if qlp_shift_needed < -(2 ** 4):
739        raise ValueError("too much negative shift needed")
740
741    qlp_max = 2 ** (qlp_precision - 1) - 1
742    qlp_min = -(2 ** (qlp_precision - 1))
743    error = 0.0
744    qlp_coeffs = []
745
746    if qlp_shift_needed >= 0:
747        for lp_coeff in lp_coefficients[order - 1]:
748            error += (lp_coeff * 2 ** qlp_shift_needed)
749            qlp_coeffs.append(min(max(int(round(error)), qlp_min), qlp_max))
750            error -= qlp_coeffs[-1]
751
752        return (qlp_coeffs, qlp_shift_needed)
753    else:
754        for lp_coeff in lp_coefficients[order - 1]:
755            error += (lp_coeff / 2 ** -qlp_shift_needed)
756            qlp_coeffs.append(min(max(int(round(error)), qlp_min), qlp_max))
757            error -= qlp_coeffs[-1]
758
759        return (qlp_coeffs, 0)
760
761
762def write_lpc_subframe(writer, options, wasted_bps, bits_per_sample,
763                       order, qlp_precision, qlp_shift_needed,
764                       qlp_coefficients, samples):
765    assert(order == len(qlp_coefficients))
766    assert(qlp_shift_needed >= 0)
767
768    # write subframe header
769    writer.build("1p 1u 5u", [1, order - 1])
770    if wasted_bps > 0:
771        writer.write(1, 1)
772        writer.unary(1, wasted_bps - 1)
773    else:
774        writer.write(1, 0)
775
776    # write warm-up samples
777    for sample in samples[0:order]:
778        writer.write_signed(bits_per_sample - wasted_bps, sample)
779
780    # write precision and shift-needed
781    writer.build("4u 5s", (qlp_precision - 1, qlp_shift_needed))
782
783    # write QLP coefficients
784    for qlp_coeff in qlp_coefficients:
785        writer.write_signed(qlp_precision, qlp_coeff)
786
787    # calculate residuals
788    residuals = []
789    coefficients = list(reversed(qlp_coefficients))
790
791    for (i, sample) in enumerate(samples[order:]):
792        residuals.append(sample - (sum([c * s for c, s in
793                                        zip(coefficients,
794                                            samples[i:i + order])]) >>
795                                   qlp_shift_needed))
796
797    # write residual block
798    encode_residuals(writer, options, order, len(samples), residuals)
799
800
801class CRC8(object):
802    TABLE = [0x00, 0x07, 0x0E, 0x09, 0x1C, 0x1B, 0x12, 0x15,
803             0x38, 0x3F, 0x36, 0x31, 0x24, 0x23, 0x2A, 0x2D,
804             0x70, 0x77, 0x7E, 0x79, 0x6C, 0x6B, 0x62, 0x65,
805             0x48, 0x4F, 0x46, 0x41, 0x54, 0x53, 0x5A, 0x5D,
806             0xE0, 0xE7, 0xEE, 0xE9, 0xFC, 0xFB, 0xF2, 0xF5,
807             0xD8, 0xDF, 0xD6, 0xD1, 0xC4, 0xC3, 0xCA, 0xCD,
808             0x90, 0x97, 0x9E, 0x99, 0x8C, 0x8B, 0x82, 0x85,
809             0xA8, 0xAF, 0xA6, 0xA1, 0xB4, 0xB3, 0xBA, 0xBD,
810             0xC7, 0xC0, 0xC9, 0xCE, 0xDB, 0xDC, 0xD5, 0xD2,
811             0xFF, 0xF8, 0xF1, 0xF6, 0xE3, 0xE4, 0xED, 0xEA,
812             0xB7, 0xB0, 0xB9, 0xBE, 0xAB, 0xAC, 0xA5, 0xA2,
813             0x8F, 0x88, 0x81, 0x86, 0x93, 0x94, 0x9D, 0x9A,
814             0x27, 0x20, 0x29, 0x2E, 0x3B, 0x3C, 0x35, 0x32,
815             0x1F, 0x18, 0x11, 0x16, 0x03, 0x04, 0x0D, 0x0A,
816             0x57, 0x50, 0x59, 0x5E, 0x4B, 0x4C, 0x45, 0x42,
817             0x6F, 0x68, 0x61, 0x66, 0x73, 0x74, 0x7D, 0x7A,
818             0x89, 0x8E, 0x87, 0x80, 0x95, 0x92, 0x9B, 0x9C,
819             0xB1, 0xB6, 0xBF, 0xB8, 0xAD, 0xAA, 0xA3, 0xA4,
820             0xF9, 0xFE, 0xF7, 0xF0, 0xE5, 0xE2, 0xEB, 0xEC,
821             0xC1, 0xC6, 0xCF, 0xC8, 0xDD, 0xDA, 0xD3, 0xD4,
822             0x69, 0x6E, 0x67, 0x60, 0x75, 0x72, 0x7B, 0x7C,
823             0x51, 0x56, 0x5F, 0x58, 0x4D, 0x4A, 0x43, 0x44,
824             0x19, 0x1E, 0x17, 0x10, 0x05, 0x02, 0x0B, 0x0C,
825             0x21, 0x26, 0x2F, 0x28, 0x3D, 0x3A, 0x33, 0x34,
826             0x4E, 0x49, 0x40, 0x47, 0x52, 0x55, 0x5C, 0x5B,
827             0x76, 0x71, 0x78, 0x7F, 0x6A, 0x6D, 0x64, 0x63,
828             0x3E, 0x39, 0x30, 0x37, 0x22, 0x25, 0x2C, 0x2B,
829             0x06, 0x01, 0x08, 0x0F, 0x1A, 0x1D, 0x14, 0x13,
830             0xAE, 0xA9, 0xA0, 0xA7, 0xB2, 0xB5, 0xBC, 0xBB,
831             0x96, 0x91, 0x98, 0x9F, 0x8A, 0x8D, 0x84, 0x83,
832             0xDE, 0xD9, 0xD0, 0xD7, 0xC2, 0xC5, 0xCC, 0xCB,
833             0xE6, 0xE1, 0xE8, 0xEF, 0xFA, 0xFD, 0xF4, 0xF3]
834
835    def __init__(self):
836        self.value = 0
837
838    def __int__(self):
839        return self.value
840
841    def update(self, byte):
842        self.value = self.TABLE[self.value ^ byte]
843
844
845class CRC16(CRC8):
846    TABLE = [0x0000, 0x8005, 0x800f, 0x000a, 0x801b, 0x001e, 0x0014, 0x8011,
847             0x8033, 0x0036, 0x003c, 0x8039, 0x0028, 0x802d, 0x8027, 0x0022,
848             0x8063, 0x0066, 0x006c, 0x8069, 0x0078, 0x807d, 0x8077, 0x0072,
849             0x0050, 0x8055, 0x805f, 0x005a, 0x804b, 0x004e, 0x0044, 0x8041,
850             0x80c3, 0x00c6, 0x00cc, 0x80c9, 0x00d8, 0x80dd, 0x80d7, 0x00d2,
851             0x00f0, 0x80f5, 0x80ff, 0x00fa, 0x80eb, 0x00ee, 0x00e4, 0x80e1,
852             0x00a0, 0x80a5, 0x80af, 0x00aa, 0x80bb, 0x00be, 0x00b4, 0x80b1,
853             0x8093, 0x0096, 0x009c, 0x8099, 0x0088, 0x808d, 0x8087, 0x0082,
854             0x8183, 0x0186, 0x018c, 0x8189, 0x0198, 0x819d, 0x8197, 0x0192,
855             0x01b0, 0x81b5, 0x81bf, 0x01ba, 0x81ab, 0x01ae, 0x01a4, 0x81a1,
856             0x01e0, 0x81e5, 0x81ef, 0x01ea, 0x81fb, 0x01fe, 0x01f4, 0x81f1,
857             0x81d3, 0x01d6, 0x01dc, 0x81d9, 0x01c8, 0x81cd, 0x81c7, 0x01c2,
858             0x0140, 0x8145, 0x814f, 0x014a, 0x815b, 0x015e, 0x0154, 0x8151,
859             0x8173, 0x0176, 0x017c, 0x8179, 0x0168, 0x816d, 0x8167, 0x0162,
860             0x8123, 0x0126, 0x012c, 0x8129, 0x0138, 0x813d, 0x8137, 0x0132,
861             0x0110, 0x8115, 0x811f, 0x011a, 0x810b, 0x010e, 0x0104, 0x8101,
862             0x8303, 0x0306, 0x030c, 0x8309, 0x0318, 0x831d, 0x8317, 0x0312,
863             0x0330, 0x8335, 0x833f, 0x033a, 0x832b, 0x032e, 0x0324, 0x8321,
864             0x0360, 0x8365, 0x836f, 0x036a, 0x837b, 0x037e, 0x0374, 0x8371,
865             0x8353, 0x0356, 0x035c, 0x8359, 0x0348, 0x834d, 0x8347, 0x0342,
866             0x03c0, 0x83c5, 0x83cf, 0x03ca, 0x83db, 0x03de, 0x03d4, 0x83d1,
867             0x83f3, 0x03f6, 0x03fc, 0x83f9, 0x03e8, 0x83ed, 0x83e7, 0x03e2,
868             0x83a3, 0x03a6, 0x03ac, 0x83a9, 0x03b8, 0x83bd, 0x83b7, 0x03b2,
869             0x0390, 0x8395, 0x839f, 0x039a, 0x838b, 0x038e, 0x0384, 0x8381,
870             0x0280, 0x8285, 0x828f, 0x028a, 0x829b, 0x029e, 0x0294, 0x8291,
871             0x82b3, 0x02b6, 0x02bc, 0x82b9, 0x02a8, 0x82ad, 0x82a7, 0x02a2,
872             0x82e3, 0x02e6, 0x02ec, 0x82e9, 0x02f8, 0x82fd, 0x82f7, 0x02f2,
873             0x02d0, 0x82d5, 0x82df, 0x02da, 0x82cb, 0x02ce, 0x02c4, 0x82c1,
874             0x8243, 0x0246, 0x024c, 0x8249, 0x0258, 0x825d, 0x8257, 0x0252,
875             0x0270, 0x8275, 0x827f, 0x027a, 0x826b, 0x026e, 0x0264, 0x8261,
876             0x0220, 0x8225, 0x822f, 0x022a, 0x823b, 0x023e, 0x0234, 0x8231,
877             0x8213, 0x0216, 0x021c, 0x8219, 0x0208, 0x820d, 0x8207, 0x0202]
878
879    def update(self, byte):
880        self.value = ((self.TABLE[(self.value >> 8) ^ byte] ^
881                       (self.value << 8)) & 0xFFFF)
882