190 lines
8.0 KiB
C
190 lines
8.0 KiB
C
////////////////////////////////////////////////////////////////////////////
|
|
// **** ADPCM-XQ **** //
|
|
// Xtreme Quality ADPCM Encoder/Decoder //
|
|
// Copyright (c) 2024 David Bryant. //
|
|
// All Rights Reserved. //
|
|
// Distributed under the BSD Software License (see license.txt) //
|
|
////////////////////////////////////////////////////////////////////////////
|
|
|
|
// adpcm-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 <stdio.h>
|
|
|
|
#include "adpcm_code.h"
|
|
|
|
#define FILTER_LENGTH 15
|
|
#define WINDOW_LENGTH 101
|
|
#define MIN_BLOCK_SAMPLES 16
|
|
|
|
static void win_average_buffer (float *samples, int sample_count, int half_width);
|
|
|
|
// Generate the shaping values for the specified buffer of stereo or mono samples,
|
|
// one shaping value output for each sample (or stereo pair of samples). This is
|
|
// calculated by filtering the audio at fs/6 (7350 Hz at 44.1 kHz) and comparing
|
|
// the averaged levels above and below that frequency. The output shaping values
|
|
// are nominally in the range of +/-1024, with 1024 indicating first-order HF boost
|
|
// shaping and -1024 for similar LF boost. However, since -1024 would result in
|
|
// infinite DC boost (not useful) a "min_value" is passed in. An output value of
|
|
// zero represents no noise shaping. For stereo input data the channels are summed
|
|
// for the calculation and the output is still just mono. Note that at the ends of
|
|
// the buffer the values diverge from true because not all the required source
|
|
// samples are visible. Use this formula to calculate the number of samples
|
|
// required for this process to "settle":
|
|
//
|
|
// int settle_distance = (WINDOW_LENGTH >> 1) + (FILTER_LENGTH >> 1) + 1;
|
|
//
|
|
// We also pass in a "last_value" so that we can smoothly interpolate from that
|
|
// to the first calculated value during the initial "unknown" samples. This
|
|
// reduces discontinuities.
|
|
|
|
void generate_dns_values (const int16_t *samples, int sample_count, int num_chans, int sample_rate,
|
|
int16_t *values, int16_t min_value, int16_t last_value)
|
|
{
|
|
float dB_offset = 7.3, dB_scaler = 64.0, max_dB, min_dB, max_ratio, min_ratio;
|
|
int filtered_count = sample_count - FILTER_LENGTH + 1, i;
|
|
float *low_freq, *high_freq;
|
|
|
|
(void) sample_rate; // unused.
|
|
memset (values, 0, sample_count * sizeof (values [0]));
|
|
|
|
if (filtered_count <= 0)
|
|
return;
|
|
|
|
low_freq = malloc (filtered_count * sizeof (float));
|
|
high_freq = malloc (filtered_count * sizeof (float));
|
|
|
|
// First, directly calculate the lowpassed audio using the 15-tap filter. This is
|
|
// a basic sinc with Hann windowing (for a fast transition) and because the filter
|
|
// is set to exactly fs/6, some terms are zero (which we can skip). Also, because
|
|
// it's linear-phase and has an odd number of terms, we can just subtract the LF
|
|
// result from the original to get the HF values.
|
|
|
|
if (num_chans == 1)
|
|
for (i = 0; i < filtered_count; ++i, ++samples) {
|
|
float filter_sum =
|
|
((int32_t) samples [0] + samples [14]) * 0.00150031 +
|
|
((int32_t) samples [2] + samples [12]) * -0.01703392 +
|
|
((int32_t) samples [3] + samples [11]) * -0.03449186 +
|
|
((int32_t) samples [5] + samples [ 9]) * 0.11776258 +
|
|
((int32_t) samples [6] + samples [ 8]) * 0.26543272 +
|
|
(int32_t) samples [7] * 0.33366033;
|
|
|
|
high_freq [i] = samples [FILTER_LENGTH >> 1] - filter_sum;
|
|
low_freq [i] = filter_sum;
|
|
}
|
|
else
|
|
for (i = 0; i < filtered_count; ++i, samples += 2) {
|
|
float filter_sum =
|
|
((int32_t) samples [ 0] + samples [ 1] + samples [28] + samples [29]) * 0.00150031 +
|
|
((int32_t) samples [ 4] + samples [ 5] + samples [24] + samples [25]) * -0.01703392 +
|
|
((int32_t) samples [ 6] + samples [ 7] + samples [22] + samples [23]) * -0.03449186 +
|
|
((int32_t) samples [10] + samples [11] + samples [18] + samples [19]) * 0.11776258 +
|
|
((int32_t) samples [12] + samples [13] + samples [16] + samples [17]) * 0.26543272 +
|
|
((int32_t) samples [14] + samples [15]) * 0.33366033;
|
|
|
|
high_freq [i] = samples [FILTER_LENGTH & ~1] + samples [FILTER_LENGTH] - filter_sum;
|
|
low_freq [i] = filter_sum;
|
|
}
|
|
|
|
// Apply a simple first-order "delta" filter to the lowpass because frequencies below fs/6
|
|
// become progressively less important for our purposes as the decorrelation filters make
|
|
// those frequencies less and less relevant. Note that after all this filtering, the
|
|
// magnitude level of the high frequency array will be 8.7 dB greater than the low frequency
|
|
// array when the filters are presented with pure white noise (determined empirically).
|
|
|
|
for (i = filtered_count - 1; i; --i)
|
|
low_freq [i] -= low_freq [i - 1];
|
|
|
|
low_freq [0] = low_freq [1]; // simply duplicate for the "unknown" sample
|
|
|
|
// Next we determine the averaged (absolute) levels for each sample using a box filter.
|
|
|
|
win_average_buffer (low_freq, filtered_count, WINDOW_LENGTH >> 1);
|
|
win_average_buffer (high_freq, filtered_count, WINDOW_LENGTH >> 1);
|
|
|
|
// calculate the minimum and maximum ratios that won't be clipped so that we only
|
|
// have to compute the logarithm when needed
|
|
|
|
max_dB = 1024 / dB_scaler - dB_offset;
|
|
min_dB = min_value / dB_scaler - dB_offset;
|
|
max_ratio = pow (10.0, max_dB / 20.0);
|
|
min_ratio = pow (10.0, min_dB / 20.0);
|
|
|
|
for (i = 0; i < filtered_count; ++i)
|
|
if (high_freq [i] > 1.0 && low_freq [i] > 1.0) {
|
|
float ratio = high_freq [i] / low_freq [i];
|
|
int shaping_value;
|
|
|
|
if (ratio >= max_ratio)
|
|
shaping_value = 1024;
|
|
else if (ratio <= min_ratio)
|
|
shaping_value = min_value;
|
|
else
|
|
shaping_value = (int) floor ((log10 (ratio) * 20.0 + dB_offset) * dB_scaler + 0.5);
|
|
|
|
values [i + (FILTER_LENGTH >> 1)] = shaping_value;
|
|
}
|
|
|
|
// interpolate the first 7 values from the supplied "last_value" to the first new value
|
|
|
|
for (i = 0; i < FILTER_LENGTH >> 1; ++i)
|
|
values [i] =
|
|
(
|
|
(int32_t) values [FILTER_LENGTH >> 1] * (i + 1) +
|
|
(int32_t) last_value * ((FILTER_LENGTH >> 1) - i) +
|
|
(FILTER_LENGTH >> 2)
|
|
) / ((FILTER_LENGTH >> 1) + 1);
|
|
|
|
// finally, copy the value at the end into the 7 final positions because unfortunately
|
|
// we have no "next_value" to interpolate with
|
|
|
|
for (i = filtered_count + (FILTER_LENGTH >> 1); i < sample_count; ++i)
|
|
values [i] = values [(FILTER_LENGTH >> 1) + filtered_count - 1];
|
|
|
|
free (low_freq);
|
|
free (high_freq);
|
|
}
|
|
|
|
// Given a buffer of floating values, apply a simple box filter of specified half width
|
|
// (total filter width is always odd) to determine the averaged magnitude at each point.
|
|
// For the ends, we use only the visible samples.
|
|
|
|
static void win_average_buffer (float *samples, int sample_count, int half_width)
|
|
{
|
|
float *output = malloc (sample_count * sizeof (float));
|
|
double sum = 0.0;
|
|
int m = 0, n = 0;
|
|
int i, j, k;
|
|
|
|
for (i = 0; i < sample_count; ++i) {
|
|
k = i + half_width + 1;
|
|
j = i - half_width;
|
|
|
|
if (k > sample_count) k = sample_count;
|
|
if (j < 0) j = 0;
|
|
|
|
while (m < j) {
|
|
if ((sum -= samples [m] * samples [m]) < 0.0) sum = 0.0;
|
|
m++;
|
|
}
|
|
|
|
while (n < k) {
|
|
sum += samples [n] * samples [n];
|
|
n++;
|
|
}
|
|
|
|
output [i] = sqrt (sum / (n - m));
|
|
}
|
|
|
|
memcpy (samples, output, sample_count * sizeof (float));
|
|
free (output);
|
|
}
|