rap 0.1.0

Audio processing in rust
/*
 *  Copyright © 2018 Gianmarco Garrisi
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

//! This crate provides audio processing features like computing the spectrogram and other stuff

extern crate audrey;
extern crate dft;
extern crate nalgebra;
extern crate num_traits;

use audrey::Reader;
use dft::{Operation, Plan, Transform};
use nalgebra::{DMatrix, DVector, Vector};

/// Compute a hamming window of the given size
fn hamming(frame_length: usize) -> DVector<f64> {
    use std::f64::consts::PI;
    let a0 = 25.0 / 46.0;
    DVector::from_iterator(
        frame_length,
        (0..frame_length)
            .map(|x| (a0 - (1.0 - a0) * ((2.0 * PI * x as f64 / (frame_length - 1) as f64).cos()))),
    )
}

/// Computes the Short Time Fourier Transform
pub fn magn_stft(
    frame_length: usize,
    frame_shift_length: usize,
    sampling_frequency: usize,
    signal: DVector<f64>,
) -> DMatrix<f64> {
    // generate a hamming window
    let window = hamming(frame_length);
    DMatrix::from_columns(
        &signal
            .as_slice()
            .windows(frame_length)
            .step_by(frame_shift_length)
            .map(|frame| {
                // multiply the frame by a hamming window
                let mut t = DVector::from_column_slice(frame_length, frame).component_mul(&window);
                // compute the fourier transform
                let plan = Plan::new(Operation::Forward, t.len());
                t.as_mut_slice().transform(&plan);
                // compute the absolute value
                t.abs()
            }).collect::<Vec<_>>(),
    )
}

/// Computes the spectrogram of the audio file.
pub fn spectrogram<T: std::io::Read + std::io::Seek>(
    mut file: Reader<T>,
    frame_length: usize,
    frame_shift_length: usize,
    sampling_frequency: usize,
) -> DMatrix<f64> {
    //TODO: find the type of Vec in a smarter way
    let v = file.samples().map(|t| t.unwrap()).collect::<Vec<f64>>();
    let signal = DVector::from_iterator(v.len(), v);
    let mut spec = magn_stft(frame_length, frame_shift_length, sampling_frequency, signal);
    spec.apply(|x| (x.min(10E-6)).log10()*10.0);
    spec
}

#[cfg(test)]
mod tests {
    use audrey::hound::WavWriter;
    use audrey::read::{open, Reader};

    use nalgebra::DVector;

    #[test]
    fn input_output() {
        let f = open("e001.wav").unwrap();
        if let Reader::Wav(mut wr) = f {
            let mut writer = WavWriter::create("e002.wav", wr.spec()).unwrap();
            for s in wr.samples::<i32>() {
                writer.write_sample(s.unwrap()).unwrap();
            }
            writer.finalize().unwrap();
        }
        assert_eq!(2 + 2, 4);
    }
}