Better Filters

All Windows calculated with pyfda (Python Filter Design Analysis Tool)
https://github.com/chipmuenk/pyfda
Window = Kaiser
beta = 8.6 (Similar to a Blackman Window)
fc = 22.5kHz
-86dB by 23kHz

This also gets rid of Linear Interpolation which leaves only Low and High both being Windowed Sinc.
This commit is contained in:
JasonLG1979 2023-07-04 04:02:39 -05:00
parent ac68d2431e
commit e31293cc10
5 changed files with 1234 additions and 246 deletions

View file

@ -1,7 +1,14 @@
use std::{mem, str::FromStr, time::Duration}; use std::{mem, str::FromStr, time::Duration};
pub use crate::dither::{mk_ditherer, DithererBuilder, TriangularDitherer}; pub use crate::dither::{mk_ditherer, DithererBuilder, TriangularDitherer};
use crate::{convert::i24, RESAMPLER_INPUT_SIZE, SAMPLE_RATE};
use crate::{
convert::i24,
filter_coefficients::{
HZ48000_HIGH, HZ48000_LOW, HZ88200_HIGH, HZ88200_LOW, HZ96000_HIGH, HZ96000_LOW,
},
RESAMPLER_INPUT_SIZE, SAMPLE_RATE,
};
// Reciprocals allow us to multiply instead of divide during interpolation. // Reciprocals allow us to multiply instead of divide during interpolation.
const HZ48000_RESAMPLE_FACTOR_RECIPROCAL: f64 = SAMPLE_RATE as f64 / 48_000.0; const HZ48000_RESAMPLE_FACTOR_RECIPROCAL: f64 = SAMPLE_RATE as f64 / 48_000.0;
@ -26,21 +33,9 @@ const HZ88200_INTERPOLATION_OUTPUT_SIZE: usize =
const HZ96000_INTERPOLATION_OUTPUT_SIZE: usize = const HZ96000_INTERPOLATION_OUTPUT_SIZE: usize =
(RESAMPLER_INPUT_SIZE as f64 * (1.0 / HZ96000_RESAMPLE_FACTOR_RECIPROCAL)) as 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;
const BLACKMAN_A2: f64 = 0.08;
// Constants for calculations
const TWO_TIMES_PI: f64 = 2.0 * std::f64::consts::PI;
const FOUR_TIMES_PI: f64 = 4.0 * std::f64::consts::PI;
#[derive(Clone, Copy, Debug, Default)] #[derive(Clone, Copy, Debug, Default)]
pub enum InterpolationQuality { pub enum InterpolationQuality {
Low, Low,
Medium,
#[default] #[default]
High, High,
} }
@ -53,7 +48,6 @@ impl FromStr for InterpolationQuality {
match s.to_lowercase().as_ref() { match s.to_lowercase().as_ref() {
"low" => Ok(Low), "low" => Ok(Low),
"medium" => Ok(Medium),
"high" => Ok(High), "high" => Ok(High),
_ => Err(()), _ => Err(()),
} }
@ -66,52 +60,24 @@ impl std::fmt::Display for InterpolationQuality {
match self { match self {
Low => write!(f, "Low"), Low => write!(f, "Low"),
Medium => write!(f, "Medium"),
High => write!(f, "High"), High => write!(f, "High"),
} }
} }
} }
impl InterpolationQuality { impl InterpolationQuality {
pub fn get_interpolation_coefficients(&self, resample_factor_reciprocal: f64) -> Vec<f64> { pub fn get_interpolation_coefficients(
let interpolation_coefficients_length = self.get_interpolation_coefficients_length(); &self,
mut coefficients: Vec<f64>,
let mut coefficients = Vec::with_capacity(interpolation_coefficients_length); resample_factor_reciprocal: f64,
) -> Vec<f64> {
if interpolation_coefficients_length == 0 {
warn!("InterpolationQuality::Low::get_interpolation_coefficients always returns an empty Vec<f64>");
warn!("Linear Interpolation does not use coefficients");
return coefficients;
}
let last_index = interpolation_coefficients_length as f64 - 1.0;
let sinc_center = last_index * 0.5;
let mut coefficient_sum = 0.0; let mut coefficient_sum = 0.0;
coefficients.extend((0..interpolation_coefficients_length).map( for (index, coefficient) in coefficients.iter_mut().enumerate() {
|interpolation_coefficient_index| { *coefficient *= Self::sinc((index as f64 * resample_factor_reciprocal).fract());
let index_float = interpolation_coefficient_index as f64;
let sample_index_fractional = (index_float * resample_factor_reciprocal).fract();
let sample_index_fractional_sinc_weight = Self::sinc(sample_index_fractional); coefficient_sum += *coefficient;
}
let fir_filter = Self::fir_filter(
index_float,
last_index,
sinc_center,
resample_factor_reciprocal,
);
let coefficient = sample_index_fractional_sinc_weight * fir_filter;
coefficient_sum += coefficient;
coefficient
},
));
coefficients coefficients
.iter_mut() .iter_mut()
@ -120,53 +86,6 @@ impl InterpolationQuality {
coefficients 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 {
Low => 0,
Medium => 129,
High => 257,
}
}
fn sinc(x: f64) -> f64 { fn sinc(x: f64) -> f64 {
if x.abs() < f64::EPSILON { if x.abs() < f64::EPSILON {
1.0 1.0
@ -175,35 +94,6 @@ impl InterpolationQuality {
pi_x.sin() / pi_x 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)] #[derive(Clone, Copy, Debug, Default)]
@ -262,10 +152,12 @@ impl std::fmt::Display for SampleRate {
} }
} }
#[derive(Clone, Copy, Debug, Default)] #[derive(Debug, Default)]
pub struct ResampleSpec { pub struct ResampleSpec {
pub resample_factor_reciprocal: f64, pub resample_factor_reciprocal: f64,
pub interpolation_output_size: usize, pub interpolation_output_size: usize,
pub high_coefficients: Vec<f64>,
pub low_coefficients: Vec<f64>,
} }
impl SampleRate { impl SampleRate {
@ -328,19 +220,27 @@ impl SampleRate {
ResampleSpec { ResampleSpec {
resample_factor_reciprocal: 1.0, resample_factor_reciprocal: 1.0,
interpolation_output_size: RESAMPLER_INPUT_SIZE, interpolation_output_size: RESAMPLER_INPUT_SIZE,
high_coefficients: vec![],
low_coefficients: vec![],
} }
} }
Hz48000 => ResampleSpec { Hz48000 => ResampleSpec {
resample_factor_reciprocal: HZ48000_RESAMPLE_FACTOR_RECIPROCAL, resample_factor_reciprocal: HZ48000_RESAMPLE_FACTOR_RECIPROCAL,
interpolation_output_size: HZ48000_INTERPOLATION_OUTPUT_SIZE, interpolation_output_size: HZ48000_INTERPOLATION_OUTPUT_SIZE,
high_coefficients: HZ48000_HIGH.to_vec(),
low_coefficients: HZ48000_LOW.to_vec(),
}, },
Hz88200 => ResampleSpec { Hz88200 => ResampleSpec {
resample_factor_reciprocal: HZ88200_RESAMPLE_FACTOR_RECIPROCAL, resample_factor_reciprocal: HZ88200_RESAMPLE_FACTOR_RECIPROCAL,
interpolation_output_size: HZ88200_INTERPOLATION_OUTPUT_SIZE, interpolation_output_size: HZ88200_INTERPOLATION_OUTPUT_SIZE,
high_coefficients: HZ88200_HIGH.to_vec(),
low_coefficients: HZ88200_LOW.to_vec(),
}, },
Hz96000 => ResampleSpec { Hz96000 => ResampleSpec {
resample_factor_reciprocal: HZ96000_RESAMPLE_FACTOR_RECIPROCAL, resample_factor_reciprocal: HZ96000_RESAMPLE_FACTOR_RECIPROCAL,
interpolation_output_size: HZ96000_INTERPOLATION_OUTPUT_SIZE, interpolation_output_size: HZ96000_INTERPOLATION_OUTPUT_SIZE,
high_coefficients: HZ96000_HIGH.to_vec(),
low_coefficients: HZ96000_LOW.to_vec(),
}, },
} }
} }

File diff suppressed because it is too large Load diff

View file

@ -13,6 +13,7 @@ pub mod config;
pub mod convert; pub mod convert;
pub mod decoder; pub mod decoder;
pub mod dither; pub mod dither;
pub mod filter_coefficients;
pub mod mixer; pub mod mixer;
pub mod normaliser; pub mod normaliser;
pub mod player; pub mod player;

View file

@ -1,6 +1,5 @@
use std::{ use std::{
collections::{vec_deque, VecDeque}, collections::{vec_deque, VecDeque},
marker::Send,
process::exit, process::exit,
sync::atomic::Ordering, sync::atomic::Ordering,
sync::mpsc, sync::mpsc,
@ -8,7 +7,7 @@ use std::{
}; };
use crate::{ use crate::{
config::{InterpolationQuality, SampleRate, NUM_FIR_FILTER_TAPS}, config::{InterpolationQuality, SampleRate},
player::PLAYER_COUNTER, player::PLAYER_COUNTER,
RESAMPLER_INPUT_SIZE, SAMPLE_RATE as SOURCE_SAMPLE_RATE, RESAMPLER_INPUT_SIZE, SAMPLE_RATE as SOURCE_SAMPLE_RATE,
}; };
@ -80,16 +79,6 @@ impl ConvolutionFilter {
} }
} }
trait MonoResampler {
fn new(sample_rate: SampleRate, interpolation_quality: InterpolationQuality) -> Self
where
Self: Sized;
fn stop(&mut self);
fn get_latency_pcm(&mut self) -> u64;
fn resample(&mut self, samples: &[f64]) -> Option<Vec<f64>>;
}
struct MonoSincResampler { struct MonoSincResampler {
interpolator: ConvolutionFilter, interpolator: ConvolutionFilter,
input_buffer: Vec<f64>, input_buffer: Vec<f64>,
@ -98,18 +87,22 @@ struct MonoSincResampler {
interpolation_output_size: usize, interpolation_output_size: usize,
} }
impl MonoResampler for MonoSincResampler { impl MonoSincResampler {
fn new(sample_rate: SampleRate, interpolation_quality: InterpolationQuality) -> Self { fn new(sample_rate: SampleRate, interpolation_quality: InterpolationQuality) -> Self {
let spec = sample_rate.get_resample_spec(); let spec = sample_rate.get_resample_spec();
let delay_line_latency = (interpolation_quality.get_interpolation_coefficients_length() let coefficients = match interpolation_quality {
as f64 InterpolationQuality::Low => spec.low_coefficients,
* spec.resample_factor_reciprocal) as u64; InterpolationQuality::High => spec.high_coefficients,
};
let delay_line_latency =
(coefficients.len() as f64 * spec.resample_factor_reciprocal) as u64;
Self { Self {
interpolator: ConvolutionFilter::new( interpolator: ConvolutionFilter::new(
interpolation_quality interpolation_quality
.get_interpolation_coefficients(spec.resample_factor_reciprocal), .get_interpolation_coefficients(coefficients, spec.resample_factor_reciprocal),
), ),
input_buffer: Vec::with_capacity(SOURCE_SAMPLE_RATE as usize), input_buffer: Vec::with_capacity(SOURCE_SAMPLE_RATE as usize),
@ -160,82 +153,6 @@ impl MonoResampler for MonoSincResampler {
} }
} }
struct MonoLinearResampler {
fir_filter: ConvolutionFilter,
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, 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: ConvolutionFilter::new(
interpolation_quality.get_fir_filter_coefficients(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.delay_line_latency
}
fn stop(&mut self) {
self.fir_filter.clear();
self.input_buffer.clear();
}
fn resample(&mut self, samples: &[f64]) -> Option<Vec<f64>> {
self.input_buffer.extend_from_slice(samples);
let num_buffer_chunks = self.input_buffer.len().saturating_div(RESAMPLER_INPUT_SIZE);
if num_buffer_chunks == 0 {
return None;
}
let input_size = num_buffer_chunks * RESAMPLER_INPUT_SIZE;
// The size of the output after interpolation.
// We have to account for the fact that to do effective linear
// interpolation we need an extra sample to be able to throw away later.
let output_size = num_buffer_chunks * self.interpolation_output_size + 1;
let mut output = Vec::with_capacity(output_size);
output.extend((0..output_size).map(|output_index| {
let sample_index = output_index as f64 * self.resample_factor_reciprocal;
let sample_index_fractional = sample_index.fract();
let sample_index = sample_index as usize;
let sample = *self.input_buffer.get(sample_index).unwrap_or(&0.0);
let next_sample = *self.input_buffer.get(sample_index + 1).unwrap_or(&0.0);
let sample_index_fractional_complementary = 1.0 - sample_index_fractional;
sample * sample_index_fractional_complementary + next_sample * sample_index_fractional
}));
// 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)
}
}
enum ResampleTask { enum ResampleTask {
Stop, Stop,
Terminate, Terminate,
@ -249,7 +166,7 @@ struct ResampleWorker {
} }
impl ResampleWorker { impl ResampleWorker {
fn new(mut resampler: impl MonoResampler + Send + 'static, name: String) -> Self { fn new(mut resampler: MonoSincResampler, name: String) -> Self {
let (task_sender, task_receiver) = mpsc::channel(); let (task_sender, task_receiver) = mpsc::channel();
let (result_sender, result_receiver) = mpsc::channel(); let (result_sender, result_receiver) = mpsc::channel();
@ -383,29 +300,12 @@ impl StereoInterleavedResampler {
let left_thread_name = format!("resampler:{player_id}:left"); let left_thread_name = format!("resampler:{player_id}:left");
let right_thread_name = format!("resampler:{player_id}:right"); let right_thread_name = format!("resampler:{player_id}:right");
match interpolation_quality { let left = MonoSincResampler::new(sample_rate, interpolation_quality);
InterpolationQuality::Low => { let right = MonoSincResampler::new(sample_rate, interpolation_quality);
debug!("Interpolation Type: Linear");
let left = MonoLinearResampler::new(sample_rate, interpolation_quality); Resampler::Worker {
let right = MonoLinearResampler::new(sample_rate, interpolation_quality); left_resampler: ResampleWorker::new(left, left_thread_name),
right_resampler: ResampleWorker::new(right, right_thread_name),
Resampler::Worker {
left_resampler: ResampleWorker::new(left, left_thread_name),
right_resampler: ResampleWorker::new(right, right_thread_name),
}
}
_ => {
debug!("Interpolation Type: Windowed Sinc");
let left = MonoSincResampler::new(sample_rate, interpolation_quality);
let right = MonoSincResampler::new(sample_rate, interpolation_quality);
Resampler::Worker {
left_resampler: ResampleWorker::new(left, left_thread_name),
right_resampler: ResampleWorker::new(right, right_thread_name),
}
}
} }
} }
}; };

View file

@ -582,7 +582,7 @@ fn get_setup() -> Setup {
).optopt( ).optopt(
"", "",
INTERPOLATION_QUALITY, INTERPOLATION_QUALITY,
"Interpolation Quality to use if Resampling {Low|Medium|High}. Defaults to High.", "Interpolation Quality to use if Resampling {Low|High}. Defaults to High.",
"QUALITY" "QUALITY"
).optopt( ).optopt(
"", "",
@ -817,7 +817,7 @@ fn get_setup() -> Setup {
INTERPOLATION_QUALITY, INTERPOLATION_QUALITY,
"", "",
interpolation_quality, interpolation_quality,
"Low, Medium, High", "Low, High",
default_value, default_value,
); );