192 lines
7.5 KiB
C
192 lines
7.5 KiB
C
////////////////////////////////////////////////////////////////////////////
|
|
// **** WAVPACK **** //
|
|
// Hybrid Lossless Wavefile Compressor //
|
|
// Copyright (c) 1998 - 2013 Conifer Software. //
|
|
// All Rights Reserved. //
|
|
// Distributed under the BSD Software License (see license.txt) //
|
|
////////////////////////////////////////////////////////////////////////////
|
|
|
|
// pack_dns.c
|
|
|
|
// This module handles the implementation of "dynamic noise shaping" which is
|
|
// designed to move the spectrum of the quantization noise introduced by lossy
|
|
// compression up or down in frequency so that it is more likely to be masked
|
|
// by the source material.
|
|
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
#include <math.h>
|
|
|
|
#include "wavpack_local.h"
|
|
|
|
static void best_floating_line (short *values, int num_values, double *initial_y, double *final_y, short *max_error);
|
|
|
|
void dynamic_noise_shaping (WavpackContext *wpc, int32_t *buffer, int shortening_allowed)
|
|
{
|
|
WavpackStream *wps = wpc->streams [wpc->current_stream];
|
|
int32_t sample_count = wps->wphdr.block_samples;
|
|
struct decorr_pass *ap = &wps->analysis_pass;
|
|
uint32_t flags = wps->wphdr.flags;
|
|
int32_t *bptr, temp, sam;
|
|
short *swptr;
|
|
int sc;
|
|
|
|
if (!wps->num_terms && sample_count > 8) {
|
|
if (flags & MONO_DATA)
|
|
for (bptr = buffer + sample_count - 3, sc = sample_count - 2; sc--;) {
|
|
sam = (3 * bptr [1] - bptr [2]) >> 1;
|
|
temp = *bptr-- - apply_weight (ap->weight_A, sam);
|
|
update_weight (ap->weight_A, 2, sam, temp);
|
|
}
|
|
else
|
|
for (bptr = buffer + (sample_count - 3) * 2 + 1, sc = sample_count - 2; sc--;) {
|
|
sam = (3 * bptr [2] - bptr [4]) >> 1;
|
|
temp = *bptr-- - apply_weight (ap->weight_B, sam);
|
|
update_weight (ap->weight_B, 2, sam, temp);
|
|
sam = (3 * bptr [2] - bptr [4]) >> 1;
|
|
temp = *bptr-- - apply_weight (ap->weight_A, sam);
|
|
update_weight (ap->weight_A, 2, sam, temp);
|
|
}
|
|
}
|
|
|
|
if (sample_count > wps->dc.shaping_samples) {
|
|
sc = sample_count - wps->dc.shaping_samples;
|
|
swptr = wps->dc.shaping_data + wps->dc.shaping_samples;
|
|
bptr = buffer + wps->dc.shaping_samples * ((flags & MONO_DATA) ? 1 : 2);
|
|
|
|
if (flags & MONO_DATA)
|
|
while (sc--) {
|
|
sam = (3 * ap->samples_A [0] - ap->samples_A [1]) >> 1;
|
|
temp = *bptr - apply_weight (ap->weight_A, sam);
|
|
update_weight (ap->weight_A, 2, sam, temp);
|
|
ap->samples_A [1] = ap->samples_A [0];
|
|
ap->samples_A [0] = *bptr++;
|
|
*swptr++ = (ap->weight_A < 256) ? 1024 : 1536 - ap->weight_A * 2;
|
|
}
|
|
else
|
|
while (sc--) {
|
|
sam = (3 * ap->samples_A [0] - ap->samples_A [1]) >> 1;
|
|
temp = *bptr - apply_weight (ap->weight_A, sam);
|
|
update_weight (ap->weight_A, 2, sam, temp);
|
|
ap->samples_A [1] = ap->samples_A [0];
|
|
ap->samples_A [0] = *bptr++;
|
|
|
|
sam = (3 * ap->samples_B [0] - ap->samples_B [1]) >> 1;
|
|
temp = *bptr - apply_weight (ap->weight_B, sam);
|
|
update_weight (ap->weight_B, 2, sam, temp);
|
|
ap->samples_B [1] = ap->samples_B [0];
|
|
ap->samples_B [0] = *bptr++;
|
|
|
|
*swptr++ = (ap->weight_A + ap->weight_B < 512) ? 1024 : 1536 - ap->weight_A - ap->weight_B;
|
|
}
|
|
|
|
wps->dc.shaping_samples = sample_count;
|
|
}
|
|
|
|
if (wpc->wvc_flag) {
|
|
int max_allowed_error = 1000000 / wpc->ave_block_samples;
|
|
short max_error, trial_max_error;
|
|
double initial_y, final_y;
|
|
|
|
if (max_allowed_error < 128)
|
|
max_allowed_error = 128;
|
|
|
|
best_floating_line (wps->dc.shaping_data, sample_count, &initial_y, &final_y, &max_error);
|
|
|
|
if (shortening_allowed && max_error > max_allowed_error) {
|
|
int min_samples = 0, max_samples = sample_count, trial_count;
|
|
double trial_initial_y, trial_final_y;
|
|
|
|
while (1) {
|
|
trial_count = (min_samples + max_samples) / 2;
|
|
|
|
best_floating_line (wps->dc.shaping_data, trial_count, &trial_initial_y,
|
|
&trial_final_y, &trial_max_error);
|
|
|
|
if (trial_max_error < max_allowed_error) {
|
|
max_error = trial_max_error;
|
|
min_samples = trial_count;
|
|
initial_y = trial_initial_y;
|
|
final_y = trial_final_y;
|
|
}
|
|
else
|
|
max_samples = trial_count;
|
|
|
|
if (min_samples > 10000 || max_samples - min_samples < 2)
|
|
break;
|
|
}
|
|
|
|
sample_count = min_samples;
|
|
}
|
|
|
|
if (initial_y < -512) initial_y = -512;
|
|
else if (initial_y > 1024) initial_y = 1024;
|
|
|
|
if (final_y < -512) final_y = -512;
|
|
else if (final_y > 1024) final_y = 1024;
|
|
#if 0
|
|
error_line ("%.2f sec, sample count = %5d, max error = %3d, range = %5d, %5d, actual = %5d, %5d",
|
|
(double) wps->sample_index / wpc->config.sample_rate, sample_count, max_error,
|
|
(int) floor (initial_y), (int) floor (final_y),
|
|
wps->dc.shaping_data [0], wps->dc.shaping_data [sample_count-1]);
|
|
#endif
|
|
if (sample_count != wps->wphdr.block_samples)
|
|
wps->wphdr.block_samples = sample_count;
|
|
|
|
if (wpc->wvc_flag) {
|
|
wps->dc.shaping_acc [0] = wps->dc.shaping_acc [1] = (int32_t) floor (initial_y * 65536.0 + 0.5);
|
|
|
|
wps->dc.shaping_delta [0] = wps->dc.shaping_delta [1] =
|
|
(int32_t) floor ((final_y - initial_y) / (sample_count - 1) * 65536.0 + 0.5);
|
|
|
|
wps->dc.shaping_array = NULL;
|
|
}
|
|
else
|
|
wps->dc.shaping_array = wps->dc.shaping_data;
|
|
}
|
|
else
|
|
wps->dc.shaping_array = wps->dc.shaping_data;
|
|
}
|
|
|
|
// Given an array of integer data (in shorts), find the linear function that most closely
|
|
// represents it (based on minimum sum of absolute errors). This is returned as the double
|
|
// precision initial & final Y values of the best-fit line. The function can also optionally
|
|
// compute and return a maximum error value (as a short). Note that the ends of the resulting
|
|
// line may fall way outside the range of input values, so some sort of clipping may be
|
|
// needed.
|
|
|
|
static void best_floating_line (short *values, int num_values, double *initial_y, double *final_y, short *max_error)
|
|
{
|
|
double left_sum = 0.0, right_sum = 0.0, center_x = (num_values - 1) / 2.0, center_y, m;
|
|
int i;
|
|
|
|
for (i = 0; i < num_values >> 1; ++i) {
|
|
right_sum += values [num_values - i - 1];
|
|
left_sum += values [i];
|
|
}
|
|
|
|
if (num_values & 1) {
|
|
right_sum += values [num_values >> 1] * 0.5;
|
|
left_sum += values [num_values >> 1] * 0.5;
|
|
}
|
|
|
|
center_y = (right_sum + left_sum) / num_values;
|
|
m = (right_sum - left_sum) / ((double) num_values * num_values) * 4.0;
|
|
|
|
if (initial_y)
|
|
*initial_y = center_y - m * center_x;
|
|
|
|
if (final_y)
|
|
*final_y = center_y + m * center_x;
|
|
|
|
if (max_error) {
|
|
double max = 0.0;
|
|
|
|
for (i = 0; i < num_values; ++i)
|
|
if (fabs (values [i] - (center_y + (i - center_x) * m)) > max)
|
|
max = fabs (values [i] - (center_y + (i - center_x) * m));
|
|
|
|
*max_error = (short) floor (max + 0.5);
|
|
}
|
|
}
|