From 3fd82de02c887b5050e1be66002e16495410dd52 Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 29 Sep 2016 22:45:59 -0700 Subject: [PATCH] ultralight LPC very nearly lossless encoder --- kitsune/codec/lpc.c | 248 ++++++++++++++++++++++++++++++++++++++ kitsune/codec/lpc.h | 17 +++ kitsune/codec/rice.c | 191 +++++++++++++++++++++++++++++ kitsune/codec/rice.h | 26 ++++ kitsune/fatfs_cmd.c | 4 +- kitsune/hlo_audio_tools.c | 142 ++++++++++++++++++++++ 6 files changed, 626 insertions(+), 2 deletions(-) create mode 100644 kitsune/codec/lpc.c create mode 100644 kitsune/codec/lpc.h create mode 100644 kitsune/codec/rice.c create mode 100644 kitsune/codec/rice.h diff --git a/kitsune/codec/lpc.c b/kitsune/codec/lpc.c new file mode 100644 index 00000000..9db40aa3 --- /dev/null +++ b/kitsune/codec/lpc.c @@ -0,0 +1,248 @@ +#include +#include +#include +#include + +#include "lpc.h" + +//Lookup Tables for 1st, 2nd and higher order quantized reflection coeffs +static const int32_t +lookup_1st_order_coeffs[128] = +{-2147483648,-2147221504,-2146435072,-2145124352,-2143289344,-2140930048,-2138046464,-2134638592,-2130706432,-2126249984,-2121269248,-2115764224,-2109734912,-2103181312,-2096103424,-2088501248,-2080374784,-2071724032,-2062548992,-2052849664,-2042626048,-2031878144,-2020605952,-2008809472,-1996488704,-1983643648,-1970274304,-1956380672,-1941962752,-1927020544,-1911554048,-1895563264,-1879048192,-1862008832,-1844445184,-1826357248,-1807745024,-1788608512,-1768947712,-1748762624,-1728053248,-1706819584,-1685061632,-1662779392,-1639972864,-1616642048,-1592786944,-1568407552,-1543503872,-1518075904,-1492123648,-1465647104,-1438646272,1411121152,-1383071744,-1354498048,-1325400064,-1295777792,-1265631232,-1234960384,-1203765248,-1172045824,-1139802112,-1107034112,-1073741824,-1039925248,-1005584384,-970719232,-935329792,-899416064,-862978048,-826015744,-788529152,-750518272,-711983104,-672923648,-633339904,-593231872,-552599552,-511442944,-469762048,-427556864,-384827392,-341573632,-297795584,-253493248,-208666624,-163315712,-117440512,-71041024,-24117248,23330816,71303168,119799808,168820736,218365952,268435456,319029248,370147328,421789696,473956352,526647296,579862528,633602048,687865856,742653952,797966336,853803008,910163968,967049216,1024458752,1082392576,1140850688,1199833088,1259339776,1319370752,1379926016,1441005568,1502609408,1564737536,1627389952,1690566656,1754267648,1818492928,1883242496,1948516352,2014314496,2080636928}; + +/*Checks if every sample of a block is equal or not*/ +int32_t check_if_constant(const int16_t *data,int32_t num_elements) +{ + int16_t temp = data[0]; + int i; + + for( i = 1; i < num_elements; i++) + { + if(temp != data[i]) + return -1; + } + + return 0; +} + +/*autocorrelation function*/ +void auto_corr_fun(int32_t *x,int32_t N,int64_t k,int16_t norm,int32_t *rxx) +{ + int64_t i, n; + int32_t sum = 0,mean = 0; + + for (i = 0; i < N; i++) + sum += x[i]; + mean = sum/N; + + for(i = 0; i <= k; i++) + { + rxx[i] = 0.0; + for (n = i; n < N; n++) + rxx[i] += ((int64_t)(x[n] - mean) * (x[n-i] - mean))>>15; + } + + if(norm) + { + for (i = 1; i <= k; i++) + rxx[i] /= rxx[0]; + rxx[0] = 1<<15; + } + + return; +} + +/*Levinson recursion algorithm*/ +void levinson(int32_t *autoc,uint8_t max_order,int32_t *ref,int32_t lpc[][MAX_LPC_ORDER]) +{ + int32_t i, j, i2; + int32_t r, err, tmp; + int32_t lpc_tmp[MAX_LPC_ORDER]; + + for(i = 0; i < max_order; i++) + lpc_tmp[i] = 0; + err = 1.0; + if(autoc) + err = autoc[0]; + + for(i = 0; i < max_order; i++) + { + if(ref) + r = ref[i]; + else + { + r = -autoc[i+1]; + for(j = 0; j < i; j++) + r -= lpc_tmp[j] * autoc[i-j]; + r /= err; + err *= 1.0 - (r * r); + } + + i2 = i >> 1; + lpc_tmp[i] = r; + for(j = 0; j < i2; j++) + { + tmp = lpc_tmp[j]; + lpc_tmp[j] += r * lpc_tmp[i-1-j]; + lpc_tmp[i-1-j] += r * tmp; + } + if(i & 1) + lpc_tmp[j] += lpc_tmp[j] * r; + + for(j = 0; j <= i; j++) + lpc[i][j] = -lpc_tmp[j]; + } + + return; +} + +/*Calculate reflection coefficients*/ +uint8_t compute_ref_coefs(int32_t *autoc,uint8_t max_order,int32_t *ref) +{ + int32_t i, j; + int32_t error; + int32_t gen[2][MAX_LPC_ORDER]; + uint8_t order_est; + + //Schurr recursion + for(i = 0; i < max_order; i++) + gen[0][i] = gen[1][i] = autoc[i+1]; + + error = autoc[0]; + ref[0] = ((int64_t)-gen[1][0]<<15) / error; + error += ((int64_t)gen[1][0] * ref[0] ) >> 15; + for(i = 1; i < max_order; i++) + { + for(j = 0; j < max_order - i; j++) + { + gen[1][j] = gen[1][j+1] + (int64_t)ref[i-1] * gen[0][j] >> 15; + gen[0][j] = (int64_t)gen[1][j+1] * ref[i-1] >> 15 + gen[0][j]; + } + ref[i] = ((int64_t)-gen[1][0]<<15) / error; + error += ((int64_t)gen[1][0] * ref[i])>>15; + } + + //Estimate optimal order using reflection coefficients + order_est = 1; + for(i = max_order - 1; i >= 0; i--) + { + if(abs(ref[i]) > 0.05*(1<<15)) + { + order_est = i+1; + break; + } + } + + return(order_est); +} + +static int32_t fastsqrt( int32_t x ) { + int64_t s = x; + x = 20000u; + x = ( x + s / x )/2; + x = ( x + s / x )/2; + x = ( x + s / x )/2; + x = ( x + s / x )/2; + return x; +} + +/*Quantize reflection coeffs*/ +int32_t qtz_ref_cof(int32_t *par,uint8_t ord,int32_t *q_ref) +{ + int i; + for( i = 0; i < ord; i++) + { + if(i == 0) + q_ref[i] = (( ((SQRT2 * fastsqrt(par[i] + 1)) >> 15) - (1<<15) )) >> 9; + } + + return(0); +} + +/*Dequantize reflection coefficients*/ +int32_t dqtz_ref_cof(const int32_t *q_ref,uint8_t ord,int32_t *ref) +{ + int i; + + if(ord <= 1) + { + ref[0] = 0; + return(0); + } + + for( i = 0; i < ord; i++) + { + if(i == 0) + ref[i] = lookup_1st_order_coeffs[q_ref[i] + 64]; + } + + return(0); +} + +/*Calculate residues from samples and lpc coefficients*/ +void calc_residue(const int32_t *samples,int64_t N,int16_t ord,int16_t Q,int64_t *coff,int32_t *residues) +{ + int64_t k, i; + int64_t corr; + int64_t y; + + corr = ((int64_t)1) << (Q - 1);//Correction term + + residues[0] = samples[0]; + for(k = 1; k <= ord; k++) + { + y = corr; + for(i = 1; i <= k; i++) + { + y += (int64_t)coff[i] * samples[k-i]; + if((k - i) == 0) + break; + } + residues[k] = samples[k] - (int32_t)(y >> Q); + } + + for(k = ord + 1; k < N; k++) + { + y = corr; + + for(i = 0; i <= ord; i++) + y += (int64_t)(coff[i] * samples[k-i]); + residues[k] = samples[k] - (int32_t)(y >> Q); + } + + return; +} + +/*Calculate samples from residues and lpc coefficients*/ +void calc_signal(const int32_t *residues,int64_t N,int16_t ord,int16_t Q,int64_t *coff,int32_t *samples) +{ + int64_t k, i; + int64_t corr; + int64_t y; + + corr = ((int64_t)1) << (Q - 1);//Correction term + + samples[0] = residues[0]; + for(k = 1; k <= ord; k++) + { + y = corr; + for(i = 1; i <= k; i++) + { + y -= (int64_t)(coff[i] * samples[k-i]); + if((k - i) == 0) + break; + } + samples[k] = residues[k] - (int32_t)(y >> Q); + } + + for(k = ord + 1; k < N; k++) + { + y = corr; + + for(i = 0; i <= ord; i++) + y -= (int64_t)(coff[i] * samples[k-i]); + samples[k] = residues[k] - (int32_t)(y >> Q); + } + + return; +} diff --git a/kitsune/codec/lpc.h b/kitsune/codec/lpc.h new file mode 100644 index 00000000..557a89a6 --- /dev/null +++ b/kitsune/codec/lpc.h @@ -0,0 +1,17 @@ +#ifndef _LPC_H_ +#define _LPC_H_ + +#define SQRT2 46340 +//1.4142135623730950488016887242096 +#define MAX_LPC_ORDER 1 + +int32_t check_if_constant(const int16_t *data,int32_t num_elements); +void auto_corr_fun(int32_t *x,int32_t N,int64_t k,int16_t norm,int32_t *rxx); +void levinson(int32_t *autoc,uint8_t max_order,int32_t *ref,int32_t lpc[][MAX_LPC_ORDER]); +uint8_t compute_ref_coefs(int32_t *autoc,uint8_t max_order,int32_t *ref); +int32_t qtz_ref_cof(int32_t *par,uint8_t ord,int32_t *q_ref); +int32_t dqtz_ref_cof(const int32_t *q_ref,uint8_t ord,int32_t *ref); +void calc_residue(const int32_t *samples,int64_t N,int16_t ord,int16_t Q,int64_t *coff,int32_t *residues); +void calc_signal(const int32_t *residues,int64_t N,int16_t ord,int16_t Q,int64_t *coff,int32_t *samples); + +#endif diff --git a/kitsune/codec/rice.c b/kitsune/codec/rice.c new file mode 100644 index 00000000..1f35143f --- /dev/null +++ b/kitsune/codec/rice.c @@ -0,0 +1,191 @@ +#include +#include +#include + +#include "rice.h" + +/* +*Golomb-Rice algorithm is optimized to work with unsigned numbers +*This function maps all positive numbers to even postive numbers +*and negative numbers to odd positive numbers +*/ +int32_t signed_to_unsigned(const uint32_t data_size,const int32_t *input,uint32_t *output) +{ + int i; + for( i = 0; i < data_size; i++) + (input[i] < 0) ? (output[i] = (-(input[i] << 1)) - 1) : (output[i] = input[i] << 1); + return(0); +} + +/* +*This function does exactly the reverse of the above function +*Maps all odd integers to negative numbers +*and even integers to positive numbers +*/ +int32_t unsigned_to_signed(const uint32_t data_size,const uint32_t *input,int32_t *output) +{ + int i; + for( i = 0; i < data_size; i++) + ((input[i] & 0x01) == 0x01) ? (output[i] = -((input[i] + 1) >> 1)) : (output[i] = input[i] >> 1); + return(0); +} + +/*calculate optimum rice parameter for encoding a block of data*/ +uint16_t get_opt_rice_param(const uint32_t *data,int32_t data_size,uint32_t *req_bits) +{ + int32_t i = 0,j; + uint32_t bits[MAX_RICE_PARAM]; + uint32_t temp; + uint16_t best_index; + + for(i = 0; i < MAX_RICE_PARAM; i++) + { + temp = 0; + for(j = 0; j < data_size; j++) + { + //data[j]/2^k; + temp += data[j] >> i; + temp += 1; + temp += i; + } + bits[i] = temp; + } + + //Get the minimum parameter + temp = bits[0]; + best_index = 0; + *req_bits = bits[0]; + + for(i = 1; i < MAX_RICE_PARAM; i++) + { + if(bits[i] < temp) + { + temp = bits[i]; + *req_bits = bits[i]; + best_index = (int16_t)i; + } + } + return best_index; +} + +/* +*Encode a block of data using Golomb-Rice coding +*For more details visit +*http://michael.dipperstein.com/rice/index.html +*/ +uint32_t rice_encode_block(int16_t param,const uint32_t *input,int32_t size,uint32_t *encoded) +{ + int32_t i,j,temp,written = 0; + uint32_t bits; + + rice_encode_context context; + rice_encode_context *ctx = &context; + + ctx->buffer = 0; + ctx->filled = 0; + ctx->bits = 0; + + for(i = 0; i < size; i++) + { + //temp = floor(input[i]/2^param) + temp = input[i] >> param; + + //Write out 'temp' number of ones + for(j = 0; j < temp; j++) + { + ctx->buffer = ctx->buffer | (0x0001 << ctx->filled); + ctx->bits++; + if(ctx->filled == 31) + { + encoded[written++] = ctx->buffer; + ctx->buffer = 0; + ctx->filled = 0; + } + else + ctx->filled++; + } + + //Write out a zero bit + ctx->buffer = ctx->buffer | (0x0000 << ctx->filled); + ctx->bits++; + if(ctx->filled == 31) + { + encoded[written++] = ctx->buffer; + ctx->buffer = 0; + ctx->filled = 0; + } + else + ctx->filled++; + + //Write out the last 'param' bits of input[i] + for (j = param - 1; j >= 0; j--) + { + ctx->buffer = ctx->buffer | (((input[i] >> j) & 0x0001) << ctx->filled); + ctx->bits++; + if(ctx->filled == 31) + { + encoded[written++] = ctx->buffer; + ctx->buffer = 0; + ctx->filled = 0; + } + else + ctx->filled++; + } + } + + //Flush buffer if not empty + if(ctx->filled != 0) + encoded[written] = ctx->buffer; + + bits = ctx->bits; + return(bits); +} + +/*Decode a block of data encoded using Golomb-Rice coding*/ +uint32_t rice_decode_block(int16_t param,const uint32_t *encoded,int32_t count,uint32_t *decoded) +{ + int32_t i = 0,j = 0,q,k = 0; + int16_t tmp; + rice_decode_context context; + rice_decode_context *ctx = &context; + ctx->filled = 0; + + while(k < count) + { + //Count ones until a zero is encountered + q = 0; + while(1) + { + if(ctx->filled == 0) + { + ctx->filled = 32; + ctx->buffer = encoded[i++]; + } + tmp = ctx->buffer & 0x0001; + ctx->buffer = ctx->buffer >> 1; + ctx->filled--; + if(tmp == 0) + break; + q++; + } + decoded[k] = q << param; //x = q * 2^k + + //Read out the last 'param' bits of input + for(j = param - 1; j >= 0; j--) + { + if(ctx->filled == 0) + { + ctx->filled = 32; + ctx->buffer = encoded[i++]; + } + tmp = ctx->buffer & 0x0001; + ctx->buffer = ctx->buffer >> 1; + ctx->filled--; + decoded[k] = decoded[k] | (tmp << j); + } + + k++; + } + + return(k); +} diff --git a/kitsune/codec/rice.h b/kitsune/codec/rice.h new file mode 100644 index 00000000..b33b19bd --- /dev/null +++ b/kitsune/codec/rice.h @@ -0,0 +1,26 @@ +#ifndef _RICE_H_ +#define _RICE_H_ + +#define MAX_RICE_PARAM 20 + +typedef struct rice_encode_context +{ + uint32_t buffer; + uint32_t filled; + uint32_t bits; +} rice_encode_context; + +typedef struct rice_decode_context +{ + uint32_t buffer; + uint32_t filled; + uint32_t bytes; +} rice_decode_context; + +int32_t signed_to_unsigned(const uint32_t data_size,const int32_t *input,uint32_t *output); +int32_t unsigned_to_signed(const uint32_t data_size,const uint32_t *input,int32_t *output); +uint16_t get_opt_rice_param(const uint32_t *data,int32_t data_size,uint32_t *req_bits); +uint32_t rice_encode_block(int16_t param,const uint32_t *input,int32_t size,uint32_t *encoded); +uint32_t rice_decode_block(int16_t param,const uint32_t *encoded,int32_t count,uint32_t *decoded); + +#endif diff --git a/kitsune/fatfs_cmd.c b/kitsune/fatfs_cmd.c index da6fea32..475e3c13 100644 --- a/kitsune/fatfs_cmd.c +++ b/kitsune/fatfs_cmd.c @@ -451,7 +451,7 @@ hlo_stream_t * open_stream_from_path(char * str, uint8_t input){ break; case 'b': case 'B': - rstr = hlo_stream_bw_limited(rstr, 32, 5000); + rstr = hlo_stream_bw_limited(rstr, 1, 5000); break; #if 0 case 'n': @@ -519,7 +519,7 @@ hlo_stream_t * open_stream_from_path(char * str, uint8_t input){ break; case 'b': case 'B': - rstr = hlo_stream_bw_limited(rstr, 32, 5000); + rstr = hlo_stream_bw_limited(rstr, 1, 5000); break; case 'i': case 'I': diff --git a/kitsune/hlo_audio_tools.c b/kitsune/hlo_audio_tools.c index b2ee622e..ddd0f54c 100644 --- a/kitsune/hlo_audio_tools.c +++ b/kitsune/hlo_audio_tools.c @@ -130,6 +130,146 @@ int hlo_filter_adpcm_encoder(hlo_stream_t * input, hlo_stream_t * output, void * } return ret; } +////------------------------------------------ +#include "codec/lpc.h" +#include "codec/rice.h" +#define BLOCK_SIZE 128 + +int hlo_filter_lossless_encoder(hlo_stream_t * input, hlo_stream_t * output, void * ctx, hlo_stream_signal signal){ + const int32_t sample_rate = 16000; + const int16_t bps = 16; + const int8_t channels = 1; + const int32_t estimated_frames = 800; //6 seconds + uint8_t opt_lpc_order = 0; + uint8_t rice_param_ref,rice_param_residue; + const char magic_number[5] = "hello"; + uint16_t req_int_ref,req_int_residues,samples_per_channel; + const int16_t Q = 35; + int32_t i,j = 0; + const uint32_t frame_sync = 0xAA55FF00; + uint32_t req_bits_ref,req_bits_residues; + const int64_t corr = ((int64_t)1) << Q; + + //Arrays + int16_t short_samples[BLOCK_SIZE]; + int32_t qtz_ref_coeffs[MAX_LPC_ORDER]; + int32_t int_samples[BLOCK_SIZE]; + int32_t residues[BLOCK_SIZE]; + uint32_t unsigned_ref[MAX_LPC_ORDER]; + uint32_t encoded_ref[MAX_LPC_ORDER]; + uint32_t u_residues[BLOCK_SIZE]; + uint32_t encoded_residues[BLOCK_SIZE]; + int64_t lpc[MAX_LPC_ORDER + 1]; + int32_t qtz_samples[BLOCK_SIZE]; + int32_t autocorr[MAX_LPC_ORDER + 1]; + int32_t ref[MAX_LPC_ORDER]; + int32_t lpc_mat[MAX_LPC_ORDER][MAX_LPC_ORDER]; + + int ret = 1; + + ret &= 0= 0){ + ret = hlo_stream_transfer_all(FROM_STREAM, input, (uint8_t*)short_samples, read_size*bps/8, 4); + if ( ret < 0 ){ + break; + } + ret /= bps/8; + + samples_per_channel = ret/channels; + + ret &= hlo_stream_transfer_all(INTO_STREAM, output, (uint8_t*)&frame_sync,sizeof(int32_t), 4); + for(i = 0; i < channels; i++) + { + //Separate channels + for(j = 0; j < samples_per_channel; j++) + short_samples[j] = buffer[channels * j + i]; + //Quantize sample data + for(j = 0; j < samples_per_channel; j++) + qtz_samples[j] = (short_samples[j]<<15); + //Calculate autocorrelation data + auto_corr_fun(qtz_samples,samples_per_channel,MAX_LPC_ORDER,1, + autocorr); + //Calculate reflection coefficients + opt_lpc_order = compute_ref_coefs(autocorr,MAX_LPC_ORDER,ref); + //Quantize reflection coefficients + qtz_ref_cof(ref,opt_lpc_order,qtz_ref_coeffs); + //signed to unsigned + signed_to_unsigned(opt_lpc_order,qtz_ref_coeffs,unsigned_ref); + + //get optimum rice param and number of bits + rice_param_ref = get_opt_rice_param(unsigned_ref,opt_lpc_order, + &req_bits_ref); + //Encode ref coeffs + req_bits_ref = rice_encode_block(rice_param_ref,unsigned_ref, + opt_lpc_order,encoded_ref); + //Determine number of ints required for storage + req_int_ref = (req_bits_ref+31)/(32); + + //Dequantize reflection + dqtz_ref_cof(qtz_ref_coeffs,opt_lpc_order,ref); + + //Reflection to lpc + levinson(NULL,opt_lpc_order,ref,lpc_mat); + + lpc[0] = 0; + for(j = 0; j < opt_lpc_order; j++) + lpc[j + 1] = corr * lpc_mat[opt_lpc_order - 1][j]; + + for(j = opt_lpc_order; j < MAX_LPC_ORDER; j++) + lpc[j + 1] = 0; + + //Copy samples + for(j = 0; j < samples_per_channel; j++) + int_samples[j] = short_samples[j]; + + //Get residues + calc_residue(int_samples,samples_per_channel,opt_lpc_order,Q,lpc, + residues); + + //signed to unsigned + signed_to_unsigned(samples_per_channel,residues,u_residues); + + //get optimum rice param and number of bits + rice_param_residue = get_opt_rice_param(u_residues, + samples_per_channel,&req_bits_residues); + + //Encode residues + req_bits_residues = rice_encode_block(rice_param_residue,u_residues, + samples_per_channel,encoded_residues); + + //Determine nunber of ints required for storage + req_int_residues = (req_bits_residues+31)/(32); + + //Write rice_params,bytes,encoded lpc_coeffs to output + ret &= 0