Audio Resample

resample.lisp

;; http://www.sony.jp/ic-recorder/sound-compare/pcm/
(defpackage resample
  (:use :common-lisp)
  (:export load-pcm
           save-pcm
           lpf
           resample))
(in-package :resample)

(defun load-pcm (path &key (signed t) (bits-per-sample 16)
                      &aux (element-type `(,(if signed 'signed-byte 'unsigned-byte) ,bits-per-sample)))
  (with-open-file (in path :element-type element-type)
    (let ((samples (make-array (file-length in) :element-type element-type)))
      (read-sequence samples in)
      samples)))


(defun save-pcm (path samples
                      &key (signed t) (bits-per-sample 16)
                      &aux (element-type `(,(if signed 'signed-byte 'unsigned-byte) ,bits-per-sample)))
  (with-open-file (out path :direction :output :if-exists :supersede :element-type element-type)
    (write-sequence samples out)
    t))

(defun resample (samples from to &key (method :simple))
  (ecase method
    (:simple (simple-resamle samples from to))
    (:upsample (upsample samples from to))
    (:linear  (linear-resample samples from to))
    (:filter (filter-resample samples from to))
    ))

(defun simple-resamle (src from to)
  (loop FOR i FROM 0 BELOW (length src) BY (/ from to)
        COLLECT (aref src (floor i))))

(defun upsample (src from to)
  (loop FOR i FROM 0 BELOW (length src) BY (/ from to)
        COLLECT
        (multiple-value-bind (x y) (floor i)
          (if (zerop y)
              (aref src x)
            0)
          )))

(defun linear-resample (src from to)
  (loop FOR i FROM 0 BELOW (length src) BY (/ from to)
        COLLECT (round (/ (+ (aref src (max 0 (1- (floor i))))
                             (aref src (floor i)))
                          2))))

;; https://github.com/lukem512/Master-the-Bass/blob/bfe278fd9f2e4a612f0b6a5018908dc58e34ac5a/trunk/src/com/masterthebass/LowPassFilter.java
(defun get-alpha (cutoff-frequency sample-rate)
  (cutoff-frequency-to-alpha cutoff-frequency sample-rate))

(defun cutoff-frequency-to-alpha (cutoff-frequency sample-rate)
  (- 1 (exp (/ (- cutoff-frequency) (/ sample-rate (* pi 2))))))

(defun lpf (samples &key (sr 96000) (cutoff 24000)) ; (prev-alpha 0))
  (let* ((ramp 0)
         (ramp-num 900)
         (alpha (get-alpha cutoff sr))
         (delta 0) ;(/ (- alpha prev-alpha) ramp-num))
         ;;(alpha prev-alpha)
         (filtered (make-array (length samples) :initial-contents samples)))
    (loop FOR i FROM 1 BELOW (length samples)
          DO
          (when (<= ramp ramp-num)
            (incf alpha delta)
            (incf ramp))
          (setf (aref filtered i)
                (+ (aref filtered (1- i))
                   (* alpha
                      (- (aref filtered i)
                         (aref filtered (1- i)))
                      )
                   ))
          )
    (map (type-of samples) (lambda (x) (round x)) filtered)))

;;http://www.ic.is.tohoku.ac.jp/~swk/lecture/yaruodsp/win.html

;; http://belle.digiweb.jp/nsystem/resize/resamp/resamp-04.html
(defun sinc (x)
  (/ (sin (* pi x)) (* pi x)))

(defparameter *dist* 30)

;; https://github.com/cbxbiker61/audiofilter/blob/33c0cc1b14a8ec5b6dc57acb2c6d01237e21a1dd/lib/dsp/Kaiser.cpp
(defun kaiser (x &key (n (* *dist* 2)) (alpha 2))
  (let ((n1 (1- n)))
    (assert (<= (abs (* x 2)) n))
    (unless (<= (abs (* x 2)) n1)
      (return-from kaiser 0))
    (/ (bessel (* pi alpha
                  (sqrt (- 1 (/ (* 4 x x) (* n1 n1))))))
       (bessel (* pi alpha)))))

;; http://www-cs.ccny.cuny.edu/~wolberg/pub/crc04.pdf
(defun kaiser-sinc (x)
  (* (sinc x) (kaiser x)))

(defun bessel (x)
  (clmath:bessel-i 0 (float x)))

(defun make-filter (scale dist p0 &key (weight-fun #'kaiser-sinc)) ;sinc))
  (let* ((step (min 1.0 scale))
         (taps (* (ceiling (/ dist step)) 2))
         (p (- (* (+ (/ taps 2) -1 p0) step))))
    (loop WITH sum = 0
          FOR i FROM 0 BELOW taps
      COLLECT
      (let ((x (cond ((= p 0)                    1)
                     ((not (<= (- dist) p dist)) 0)
                     (t                          (funcall weight-fun p)))))
        (incf sum x)
        (incf p step)
        x)
      INTO coef
      FINALLY
      (return (map 'list (lambda (x) (/ x sum)) coef)))))

(defun filter-resample (src from to)
  (let* ((dist *dist*)
         (scale (/ to from))
         (ox (/ (- (/ 1 scale) 1) 2))
         (coef-map (make-hash-table)))
    (loop FOR i FROM 0 BELOW (ceiling (* (length src) scale))
      COLLECT
      (let* ((p0 (abs (+ (/ i scale) ox)))
             (ip (floor p0))
             (coef (if #1=(gethash (- p0 ip) coef-map)
                     #1#
                     (setf #1# (make-filter scale dist (- p0 ip)))))
             (coef-len (length coef))
             (out 0))
        (loop FOR j FROM 0
              FOR c IN coef
          DO
          (let* ((index (abs (+ ip (- (/ coef-len 2)) j 1)))
                 (index (if (>= index (length src))
                            (- (* (length src) 2) index 1)
                          index)))
            (incf out (* (aref src index) c))))
        out)
      INTO out-list
      FINALLY (return (map `(simple-array ,(array-element-type src) (,(length out-list)))
                           (lambda (x)
                             (max (min (round x) (1- (ash 1 15))) (- (ash 1 15)))) ; xxx: hard-coding
                           out-list))
      )))

test.lisp

;; memo: spek: a spectrum analyzer
;; memo: サンプル音源: http://www.sony.jp/ic-recorder/sound-compare/pcm/
;; memo: サンプル音源: http://creolestream.com/359

(push #P"./" asdf:*central-registry*)
(push #P"clmath/" asdf:*central-registry*)


(subseq *src* 110000)

(subseq *dst* 110000)

(reduce #'max (map 'list (lambda (x y) (- x y))
                   (subseq *src* 0 100000)
                   *dst*
                   ))

(defparameter *bits* 16)

(defparameter *src* (resample:load-pcm "/home/ohta/pcm/mono_aroj_16.raw" :bits-per-sample *bits*))

(resample:save-pcm "/tmp/dst-kaiser.raw" (resample:resample *src* 44100 11025 :method :filter) :bits-per-sample *bits*)
(resample:save-pcm "/tmp/simple.raw" (resample:resample *src* 44100 11025 :method :simple) :bits-per-sample *bits*)

(resample:save-pcm "/tmp/dst-kaiser.raw" (resample:resample *src* 44100 48000 :method :filter) :bits-per-sample *bits*)
(resample:save-pcm "/tmp/dst-simple.raw" (resample:resample *src* 44100 48000 :method :simple) :bits-per-sample *bits*)

(resample:save-pcm "/tmp/lpf.raw"
                   (resample:resample
                    (loop with s = *src*
                          repeat 20
                          DO (setf s (resample:lpf s :sr 44100 :cutoff (floor 44100 8)))
                          FINALLY (return s))
                    44100
                    11025)
                   :bits-per-sample *bits*)

(resample:save-pcm "/tmp/dst.raw"
                   (loop with s = *src*
                         repeat 100
                         DO (setf s (resample:lpf s :sr 44100 :cutoff 11000))
                         FINALLY (return s))
                   :bits-per-sample 16)

(resample:save-pcm "/tmp/dst.raw"
                   (loop with s = (resample:resample *src* 44100 88200 :method :upsample)
                         repeat 10
                         DO (setf s (resample:lpf s :sr 88200 :cutoff 22050))
                         FINALLY (return s))
                   :bits-per-sample 16)

;; polyphase
;; http://www.ws.binghamton.edu/fowler/fowler%20personal%20page/EE521_files/IV-05%20Polyphase%20FIlters_2007.pdf
;; http://www.ee.ic.ac.uk/hp/staff/dmb/courses/DSPDF/01200_Polyphase.pdf