Anti-alias Linear Interpolation also.

This commit is contained in:
JasonLG1979 2023-07-03 14:10:21 -05:00
parent cb8f6c954d
commit 33e2ce65d3
2 changed files with 150 additions and 32 deletions

View file

@ -26,6 +26,8 @@ const HZ88200_INTERPOLATION_OUTPUT_SIZE: usize =
const HZ96000_INTERPOLATION_OUTPUT_SIZE: usize =
(RESAMPLER_INPUT_SIZE as f64 * (1.0 / HZ96000_RESAMPLE_FACTOR_RECIPROCAL)) as usize;
pub const NUM_FIR_FILTER_TAPS: usize = 5;
// Blackman Window coefficients
const BLACKMAN_A0: f64 = 0.42;
const BLACKMAN_A1: f64 = 0.5;
@ -77,7 +79,9 @@ impl InterpolationQuality {
let mut coefficients = Vec::with_capacity(interpolation_coefficients_length);
if interpolation_coefficients_length == 0 {
warn!("InterpolationQuality::Low::get_interpolation_coefficients always returns an empty Vec<f64> because Linear Interpolation does not use coefficients");
warn!("InterpolationQuality::Low::get_interpolation_coefficients always returns an empty Vec<f64>");
warn!("Linear Interpolation does not use coefficients");
return coefficients;
}
@ -91,30 +95,17 @@ impl InterpolationQuality {
|interpolation_coefficient_index| {
let index_float = interpolation_coefficient_index as f64;
let sample_index_fractional = (index_float * resample_factor_reciprocal).fract();
let sinc_center_offset =
// The resample_factor_reciprocal also happens to be our
// anti-alias cutoff. In this case it represents the minimum
// output bandwidth needed to fully represent the input.
(index_float - sinc_center) * resample_factor_reciprocal;
let sample_index_fractional_sinc_weight = Self::sinc(sample_index_fractional);
let sinc_value = Self::sinc(sinc_center_offset);
// Calculate the Blackman window function for the given center offset
// w(n) = A0 - A1*cos(2πn / (N-1)) + A2*cos(4πn / (N-1)),
// where n is the center offset, N is the window size,
// and A0, A1, A2 are precalculated coefficients
let fir_filter = Self::fir_filter(
index_float,
last_index,
sinc_center,
resample_factor_reciprocal,
);
let two_pi_n = TWO_TIMES_PI * index_float;
let four_pi_n = FOUR_TIMES_PI * index_float;
let blackman_window_value = BLACKMAN_A0
- BLACKMAN_A1 * (two_pi_n / last_index).cos()
+ BLACKMAN_A2 * (four_pi_n / last_index).cos();
let sinc_window = sinc_value * blackman_window_value;
let coefficient = sinc_window * sample_index_fractional_sinc_weight;
let coefficient = sample_index_fractional_sinc_weight * fir_filter;
coefficient_sum += coefficient;
@ -129,6 +120,44 @@ impl InterpolationQuality {
coefficients
}
pub fn get_fir_filter_coefficients(&self, resample_factor_reciprocal: f64) -> Vec<f64> {
let mut coefficients = Vec::with_capacity(NUM_FIR_FILTER_TAPS);
if self.get_interpolation_coefficients_length() != 0 {
warn!("InterpolationQuality::Medium/High::get_fir_filter_coefficients always returns an empty Vec<f64>");
warn!("The FIR Filter coefficients are a part of the Windowed Sinc Interpolation coefficients");
return coefficients;
}
let last_index = NUM_FIR_FILTER_TAPS as f64 - 1.0;
let sinc_center = last_index * 0.5;
let mut coefficient_sum = 0.0;
coefficients.extend(
(0..NUM_FIR_FILTER_TAPS).map(|fir_filter_coefficient_index| {
let coefficient = Self::fir_filter(
fir_filter_coefficient_index as f64,
last_index,
sinc_center,
resample_factor_reciprocal,
);
coefficient_sum += coefficient;
coefficient
}),
);
coefficients
.iter_mut()
.for_each(|coefficient| *coefficient /= coefficient_sum);
coefficients
}
pub fn get_interpolation_coefficients_length(&self) -> usize {
use InterpolationQuality::*;
match self {
@ -146,6 +175,35 @@ impl InterpolationQuality {
pi_x.sin() / pi_x
}
}
fn blackman(index: f64, last_index: f64) -> f64 {
// Calculate the Blackman window function for the given center offset
// w(n) = A0 - A1*cos(2πn / (N-1)) + A2*cos(4πn / (N-1)),
// where n is the center offset, N is the window size,
// and A0, A1, A2 are precalculated coefficients
let two_pi_n = TWO_TIMES_PI * index;
let four_pi_n = FOUR_TIMES_PI * index;
BLACKMAN_A0 - BLACKMAN_A1 * (two_pi_n / last_index).cos()
+ BLACKMAN_A2 * (four_pi_n / last_index).cos()
}
fn fir_filter(
index: f64,
last_index: f64,
sinc_center: f64,
resample_factor_reciprocal: f64,
) -> f64 {
// The resample_factor_reciprocal also happens to be our
// anti-alias cutoff. In this case it represents the minimum
// output bandwidth needed to fully represent the input.
let adjusted_sinc_center_offset = (index - sinc_center) * resample_factor_reciprocal;
let sinc_value = Self::sinc(adjusted_sinc_center_offset);
let blackman_window_value = Self::blackman(index, last_index);
sinc_value * blackman_window_value
}
}
#[derive(Clone, Copy, Debug, Default)]

View file

@ -8,28 +8,28 @@ use std::{
};
use crate::{
config::{InterpolationQuality, SampleRate},
config::{InterpolationQuality, SampleRate, NUM_FIR_FILTER_TAPS},
player::PLAYER_COUNTER,
RESAMPLER_INPUT_SIZE, SAMPLE_RATE as SOURCE_SAMPLE_RATE,
};
struct DelayLine {
buffer: VecDeque<f64>,
interpolation_coefficients_length: usize,
coefficients_length: usize,
}
impl DelayLine {
fn new(interpolation_coefficients_length: usize) -> DelayLine {
fn new(coefficients_length: usize) -> DelayLine {
Self {
buffer: VecDeque::with_capacity(interpolation_coefficients_length),
interpolation_coefficients_length,
buffer: VecDeque::with_capacity(coefficients_length),
coefficients_length,
}
}
fn push(&mut self, sample: f64) {
self.buffer.push_back(sample);
while self.buffer.len() > self.interpolation_coefficients_length {
while self.buffer.len() > self.coefficients_length {
self.buffer.pop_front();
}
}
@ -48,12 +48,21 @@ impl<'a> IntoIterator for &'a DelayLine {
}
}
trait Convoluter {
fn new(interpolation_quality: InterpolationQuality, resample_factor_reciprocal: f64) -> Self
where
Self: Sized;
fn convolute(&mut self, sample: f64) -> f64;
fn clear(&mut self);
}
struct WindowedSincInterpolator {
interpolation_coefficients: Vec<f64>,
delay_line: DelayLine,
}
impl WindowedSincInterpolator {
impl Convoluter for WindowedSincInterpolator {
fn new(interpolation_quality: InterpolationQuality, resample_factor_reciprocal: f64) -> Self {
let interpolation_coefficients =
interpolation_quality.get_interpolation_coefficients(resample_factor_reciprocal);
@ -66,7 +75,7 @@ impl WindowedSincInterpolator {
}
}
fn interpolate(&mut self, sample: f64) -> f64 {
fn convolute(&mut self, sample: f64) -> f64 {
// Since our interpolation coefficients are pre-calculated
// we can basically pretend like the Interpolator is a FIR filter.
self.delay_line.push(sample);
@ -85,6 +94,41 @@ impl WindowedSincInterpolator {
}
}
struct LinearFirFilter {
fir_filter_coefficients: Vec<f64>,
delay_line: DelayLine,
}
impl Convoluter for LinearFirFilter {
fn new(interpolation_quality: InterpolationQuality, resample_factor_reciprocal: f64) -> Self {
let fir_filter_coefficients =
interpolation_quality.get_fir_filter_coefficients(resample_factor_reciprocal);
let delay_line = DelayLine::new(fir_filter_coefficients.len());
Self {
fir_filter_coefficients,
delay_line,
}
}
fn convolute(&mut self, sample: f64) -> f64 {
self.delay_line.push(sample);
// Temporal convolution
self.fir_filter_coefficients
.iter()
.zip(&self.delay_line)
.fold(0.0, |acc, (coefficient, delay_line_sample)| {
acc + coefficient * delay_line_sample
})
}
fn clear(&mut self) {
self.delay_line.clear();
}
}
trait MonoResampler {
fn new(sample_rate: SampleRate, interpolation_quality: InterpolationQuality) -> Self
where
@ -156,7 +200,7 @@ impl MonoResampler for MonoSincResampler {
// out the other side is Sinc Windowed Interpolated samples.
let sample_index = (ouput_index as f64 * self.resample_factor_reciprocal) as usize;
let sample = self.input_buffer[sample_index];
self.interpolator.interpolate(sample)
self.interpolator.convolute(sample)
}));
self.input_buffer.drain(..input_size);
@ -166,27 +210,39 @@ impl MonoResampler for MonoSincResampler {
}
struct MonoLinearResampler {
fir_filter: LinearFirFilter,
input_buffer: Vec<f64>,
resample_factor_reciprocal: f64,
delay_line_latency: u64,
interpolation_output_size: usize,
}
impl MonoResampler for MonoLinearResampler {
fn new(sample_rate: SampleRate, _: InterpolationQuality) -> Self {
fn new(sample_rate: SampleRate, interpolation_quality: InterpolationQuality) -> Self {
let spec = sample_rate.get_resample_spec();
let delay_line_latency =
(NUM_FIR_FILTER_TAPS as f64 * spec.resample_factor_reciprocal) as u64;
Self {
fir_filter: LinearFirFilter::new(
interpolation_quality,
spec.resample_factor_reciprocal,
),
input_buffer: Vec::with_capacity(SOURCE_SAMPLE_RATE as usize),
resample_factor_reciprocal: spec.resample_factor_reciprocal,
delay_line_latency,
interpolation_output_size: spec.interpolation_output_size,
}
}
fn get_latency_pcm(&mut self) -> u64 {
self.input_buffer.len() as u64
self.input_buffer.len() as u64 + self.delay_line_latency
}
fn stop(&mut self) {
self.fir_filter.clear();
self.input_buffer.clear();
}
@ -220,6 +276,10 @@ impl MonoResampler for MonoLinearResampler {
// Remove the last garbage sample.
output.pop();
output
.iter_mut()
.for_each(|sample| *sample = self.fir_filter.convolute(*sample));
self.input_buffer.drain(..input_size);
Some(output)