public domain pull-based linear/hermite/sinc signal resampling code in rust

resample.rs

// A public domain rust implementation of three common signal resampling algorithms.

/*
 Released under the following license (CC0).

 Creative Commons Legal Code
 
 CC0 1.0 Universal
 Statement of Purpose
 
 The laws of most jurisdictions throughout the world automatically confer
 exclusive Copyright and Related Rights (defined below) upon the creator
 and subsequent owner(s) (each and all, an "owner") of an original work of
 authorship and/or a database (each, a "Work").
 
 Certain owners wish to permanently relinquish those rights to a Work for
 the purpose of contributing to a commons of creative, cultural and
 scientific works ("Commons") that the public can reliably and without fear
 of later claims of infringement build upon, modify, incorporate in other
 works, reuse and redistribute as freely as possible in any form whatsoever
 and for any purposes, including without limitation commercial purposes.
 These owners may contribute to the Commons to promote the ideal of a free
 culture and the further production of creative, cultural and scientific
 works, or to gain reputation or greater distribution for their Work in
 part through the use and efforts of others.
 
 For these and/or other purposes and motivations, and without any
 expectation of additional consideration or compensation, the person
 associating CC0 with a Work (the "Affirmer"), to the extent that he or she
 is an owner of Copyright and Related Rights in the Work, voluntarily
 elects to apply CC0 to the Work and publicly distribute the Work under its
 terms, with knowledge of his or her Copyright and Related Rights in the
 Work and the meaning and intended legal effect of CC0 on those rights.
 
 1. Copyright and Related Rights. A Work made available under CC0 may be
 protected by copyright and related or neighboring rights ("Copyright and
 Related Rights"). Copyright and Related Rights include, but are not
 limited to, the following:
 
   i. the right to reproduce, adapt, distribute, perform, display,
      communicate, and translate a Work;
  ii. moral rights retained by the original author(s) and/or performer(s);
 iii. publicity and privacy rights pertaining to a person's image or
      likeness depicted in a Work;
  iv. rights protecting against unfair competition in regards to a Work,
      subject to the limitations in paragraph 4(a), below;
   v. rights protecting the extraction, dissemination, use and reuse of data
      in a Work;
  vi. database rights (such as those arising under Directive 96/9/EC of the
      European Parliament and of the Council of 11 March 1996 on the legal
      protection of databases, and under any national implementation
      thereof, including any amended or successor version of such
      directive); and
 vii. other similar, equivalent or corresponding rights throughout the
      world based on applicable law or treaty, and any national
      implementations thereof.
 
 2. Waiver. To the greatest extent permitted by, but not in contravention
 of, applicable law, Affirmer hereby overtly, fully, permanently,
 irrevocably and unconditionally waives, abandons, and surrenders all of
 Affirmer's Copyright and Related Rights and associated claims and causes
 of action, whether now known or unknown (including existing as well as
 future claims and causes of action), in the Work (i) in all territories
 worldwide, (ii) for the maximum duration provided by applicable law or
 treaty (including future time extensions), (iii) in any current or future
 medium and for any number of copies, and (iv) for any purpose whatsoever,
 including without limitation commercial, advertising or promotional
 purposes (the "Waiver"). Affirmer makes the Waiver for the benefit of each
 member of the public at large and to the detriment of Affirmer's heirs and
 successors, fully intending that such Waiver shall not be subject to
 revocation, rescission, cancellation, termination, or any other legal or
 equitable action to disrupt the quiet enjoyment of the Work by the public
 as contemplated by Affirmer's express Statement of Purpose.
 
 3. Public License Fallback. Should any part of the Waiver for any reason
 be judged legally invalid or ineffective under applicable law, then the
 Waiver shall be preserved to the maximum extent permitted taking into
 account Affirmer's express Statement of Purpose. In addition, to the
 extent the Waiver is so judged Affirmer hereby grants to each affected
 person a royalty-free, non transferable, non sublicensable, non exclusive,
 irrevocable and unconditional license to exercise Affirmer's Copyright and
 Related Rights in the Work (i) in all territories worldwide, (ii) for the
 maximum duration provided by applicable law or treaty (including future
 time extensions), (iii) in any current or future medium and for any number
 of copies, and (iv) for any purpose whatsoever, including without
 limitation commercial, advertising or promotional purposes (the
 "License"). The License shall be deemed effective as of the date CC0 was
 applied by Affirmer to the Work. Should any part of the License for any
 reason be judged legally invalid or ineffective under applicable law, such
 partial invalidity or ineffectiveness shall not invalidate the remainder
 of the License, and in such case Affirmer hereby affirms that he or she
 will not (i) exercise any of his or her remaining Copyright and Related
 Rights in the Work or (ii) assert any associated claims and causes of
 action with respect to the Work, in either case contrary to Affirmer's
 express Statement of Purpose.
 
 4. Limitations and Disclaimers.
 
  a. No trademark or patent rights held by Affirmer are waived, abandoned,
     surrendered, licensed or otherwise affected by this document.
  b. Affirmer offers the Work as-is and makes no representations or
     warranties of any kind concerning the Work, express, implied,
     statutory or otherwise, including without limitation warranties of
     title, merchantability, fitness for a particular purpose, non
     infringement, or the absence of latent or other defects, accuracy, or
     the present or absence of errors, whether or not discoverable, all to
     the greatest extent permissible under applicable law.
  c. Affirmer disclaims responsibility for clearing rights of other persons
     that may apply to the Work or any use thereof, including without
     limitation any person's Copyright and Related Rights in the Work.
     Further, Affirmer disclaims responsibility for obtaining any necessary
     consents, permissions or other rights required for any use of the
     Work.
  d. Affirmer understands and acknowledges that Creative Commons is not a
     party to this document and has no duty or obligation with respect to
     this CC0 or use of the Work.
*/

// get sample from slice, falling back to 0.0 if outside range
// location is floored
fn pull(data : &[f32], sample : f32) -> f32
{
    let i = sample.floor() as isize;
    if i < 0 || i as usize >= data.len()
    {
        return 0.0;
    }
    data[i as usize]
}

// normalized cardinal sine (sinc) function
fn sinc(mut x : f32) -> f32
{
    if x == 0.0
    {
        return 1.0;
    }
    x *= std::f32::consts::PI;
    x.sin()/x
}

// Sinc interpolation. See internal comments for information about window size.
// Uses a squared welch window `(1-x^2)^2` with a radius of 16 (total size of 32).
// Ratio above 1.0 when upsampling, below 1.0 when downsampling.
fn sinc_pull(data : &[f32], sample : f32, mut ratio : f32) -> f32
{
    // Window radius (i.e. will sample approximately width*2 output samples worth of input samples)
    // Lower values give progressively more artifacts but are linearly cheaper.
    // Values less than 4.0 are pointless, because goodhermite gives better results.
    // Values less than 1.0 are invalid and will probably result in serious, harsh distortion.
    // Values less than 0.5 are DEFINITELY invalid and WILL result in REALLY bad amplitude modulation effects near source nyquist.
    let mut width = 16.0;
    // if downsampling, increase window size to do lowpass, with kernel in terms of destination samples
    // otherwise suppress ratio calculation, so that the kernel is in terms of source samples
    if ratio < 1.0
    {
        width /= ratio;
    }
    else
    {
        ratio = 1.0;
    }
    let t = sample-sample.floor();
    let mut out = 0.0;
    let mut normalize = 0.0;
    // iterate over source samples
    for _i in (-width+1.0).ceil() as isize..=(width).floor() as isize
    {
        // location relative to destination sample in terms of source samples
        let i = _i as f32 - t;
        let iw = i/width;
        // location within window in range -1 to 1
        let window = (1.0-iw*iw).powi(2);
        // calculate kernel, then window it
        let kernel = sinc(i * ratio) * window;
        normalize += kernel;
        // grab and multiply source sample, add to sum
        out += pull(data, _i as f32 + sample.floor())*kernel;
    }
    out /= normalize; // ensure that DC offset signals don't introduce higher frequency tones
    out
}
// Sinc interpolation.
// Uses a squared welch window `(1-x^2)^2`.
// With window size customization. Window size is in terms of destination samples when downsampling, source samples when upsampling.
// Ratio above 1.0 when upsampling, below 1.0 when downsampling.
fn sinc_pull_width(data : &[f32], sample : f32, mut ratio : f32, mut width : f32) -> f32
{
    // if downsampling, increase window size to do lowpass, with kernel in terms of destination samples
    // otherwise suppress ratio calculation, so that the kernel is in terms of source samples
    if ratio < 1.0
    {
        width /= ratio;
    }
    else
    {
        ratio = 1.0;
    }
    let t = sample-sample.floor();
    let mut out = 0.0;
    let mut normalize = 0.0;
    // iterate over source samples
    for _i in (-width+1.0).ceil() as isize..=(width).floor() as isize
    {
        // location relative to destination sample in terms of source samples
        let i = _i as f32 - t;
        // location within window in range -1 to 1
        let iw = i/width;
        let window = (1.0-iw*iw).powi(2);
        // calculate kernel, then window it
        let kernel = sinc(i * ratio) * window;
        normalize += kernel;
        // grab and multiply source sample, add to sum
        out += pull(data, sample.floor() + _i as f32) * kernel;
    }
    out /= normalize; // ensure that DC offset signals don't introduce higher frequency tones
    out
}

// Handles both upsampling and downsampling properly.
// Uses the cubic hermite interpolation of a simple 1.0 sample as a kernel.
// Same result as normal cubic interpolation when upsampling.
// Stretches the kernel to function as a lowpass when downsampling.
// This means that it samples more than just four-ish input samples when downsampling.
fn goodhermite_pull(data : &[f32], sample : f32, mut ratio : f32) -> f32
{
    // do not change (raising from 2.0 has no effect, lowering from 2.0 breaks the interpolation)
    let mut width = 2.0;
    // if downsampling, increase window size to do lowpass, with kernel in terms of destination samples
    // otherwise suppress ratio calculation, so that the kernel is in terms of source samples
    if ratio < 1.0
    {
        width /= ratio;
    }
    else
    {
        ratio = 1.0;
    }
    let t = sample-sample.floor();
    let mut out = 0.0;
    let mut normalize = 0.0;
    // iterate over source samples
    for _i in (-width+1.0).ceil() as isize..=(width).floor() as isize
    {
        // location relative to destination sample in terms of source samples
        let i = _i as f32 - t;
        // calculate kernel
        let kernel = hermite_kernel(i * ratio);
        normalize += kernel;
        // grab and multiply source sample, add to sum
        out += pull(data, _i as f32 + sample.floor())*kernel;
    }
    out /= normalize; // ensure that DC offset signals don't introduce higher frequency tones
    out
}
const _hermite_kernel_data : [f32; 1] = [1.0];
fn hermite_kernel(sample : f32) -> f32
{
    return hermite_pull(&_hermite_kernel_data, sample);
}

// Simple cubic interpolation (hermite spline with Catmull-Rom tangent calculation).
// Handles upsampling properly, but not downsampling. Uses Catmull-Rom tangent calculation.
fn hermite_pull(data : &[f32], sample : f32) -> f32
{
    let t = sample-sample.floor();
    let a = pull(data, sample.floor()-1.0);
    let b = pull(data, sample.floor()-0.0);
    let c = pull(data, sample.floor()+1.0);
    let d = pull(data, sample.floor()+2.0);
    
    let ta = (c-a)/2.0;
    let tb = (d-b)/2.0;
    
      b  * ( 2.0*t*t*t - 3.0*t*t + 1.0)
    + c  * (-2.0*t*t*t + 3.0*t*t)
    + ta * (     t*t*t - 2.0*t*t + t)
    + tb * (     t*t*t -     t*t)
}

// Linear interpolation.
fn linear_pull(data : &[f32], sample : f32) -> f32
{
    let t = sample-sample.floor();
    let sample_a = pull(data, sample-0.0);
    let sample_b = pull(data, sample+1.0);
    return sample_a * (1.0-t) + sample_b * t;
}

// demo function
fn main()
{
    /*
    for i in -10..=10
    {
        let a = badsinc_pull(&_hermite_kernel_data, i as f32/5.0, 5.0);
        let b = goodhermite_pull(&_hermite_kernel_data, i as f32/5.0, 5.0);
        let c = (a-b).abs();
        println!("{a} {b} {c}");
    }
    return;
    */
    
    use std::fs::File;
    use std::path::Path;
    
    let args : Vec<String> = std::env::args().collect();
    
    let mut in_file = File::open(Path::new(&args[1])).unwrap();
    let (mut header, data) = wav::read(&mut in_file).unwrap();
    
    let data = data.as_sixteen().unwrap().iter().map(|x| *x as f32 / 32768.0).collect::<Vec<_>>();
    
    header.audio_format = wav::header::WAV_FORMAT_IEEE_FLOAT;
    header.bits_per_sample = 32;
    
    let mut newdata = Vec::new();
    
    let rate_conversion = 44100.0/header.sampling_rate as f32;
    println!("{rate_conversion}");
    header.sampling_rate = 44100;
    let new_sample_count = (data.len() as f32 * rate_conversion) as usize;
    
    for i in 0..new_sample_count
    {
        let source_i = i as f32/rate_conversion;
        newdata.push(sinc_pull(&data, source_i, rate_conversion));
        //newdata.push(hermite_pull(&data, source_i));
    }
    
    let newdata = wav::bit_depth::BitDepth::ThirtyTwoFloat(newdata);
    let mut out_file = File::create(Path::new("_test_out.wav")).unwrap();
    wav::write(header, &newdata, &mut out_file).unwrap();
}