#lang typed/racket
(require/typed "rsound.rkt"
[default-sample-rate (-> Real)])
(provide response/raw
response/mag
roots->poly
coefficients->poly
poles&zeros->poly
coefficient-sets->poly
roots->coefficients
real-part/ck
lpf-coefficients
up-to-power-of-two
all-but-n
product-of
sum-of)
(: i Complex)
(define i (exact->inexact (sqrt -1)))
(define twopi (* 2 pi))
(define-type Poly (Complex -> Complex))
(define-type Z-Plane-Points (Listof Complex))
(define-type Zeros Z-Plane-Points)
(define-type Poles Z-Plane-Points)
(define-type Coefficients (Listof Real))
(define-type Frequency Nonnegative-Real)
(define-type Signal (Nonnegative-Fixnum -> Real))
(: response/raw (Poly -> Frequency -> Complex))
(define (response/raw poly)
(define sr-inv (exact->inexact
(/ 1 (default-sample-rate))))
(lambda (omega)
(define z (exp (* i sr-inv twopi omega)))
(poly z)))
(: response/mag (Poly -> Frequency -> Real))
(define ((response/mag poly) frequency)
(define mag (exact->inexact
(magnitude ((response/raw poly) frequency))))
(max -100 (* 10 (/ (log (expt mag 2)) (log 10))))
(* 10 (/ (ann (log (ann (max 1.0e-6 mag)
Positive-Real))
Real)
(/ (log 10) 2))))
(: roots->poly (Zeros -> Poly))
(define (roots->poly roots)
(coefficients->poly (roots->coefficients roots)))
(: coefficients->poly (Coefficients -> Poly))
(define ((coefficients->poly coefficients) x)
(for/fold: ([so-far : Complex 0.0+0.0i])
([coefficient (in-list coefficients)])
(+ coefficient (* so-far x))))
(: poles&zeros->poly (Poles Zeros -> Poly))
(define (poles&zeros->poly poles zeros)
(coefficient-sets->poly (roots->coefficients poles)
(roots->coefficients zeros)))
(: coefficient-sets->poly (Coefficients Coefficients -> Poly))
(define (coefficient-sets->poly fb-coefficients ff-coefficients)
(let ([feedback-poly (coefficients->poly fb-coefficients)]
[feedforward-poly (coefficients->poly ff-coefficients)])
(lambda (x)
(/ (feedforward-poly x)
(feedback-poly x)))))
(: roots->coefficients (Z-Plane-Points -> Coefficients))
(define (roots->coefficients z-plane-points)
(let ([neg-points (map - z-plane-points)])
(reverse
(for/list: ([exponent : Exact-Nonnegative-Integer
(ann
(in-range (ann (add1 (length neg-points))
Nonnegative-Fixnum))
(Sequenceof Nonnegative-Fixnum))])
(real-part/ck
(sum-of (map product-of (all-but-n exponent neg-points))))))))
(: lpf-coefficients (Real -> Coefficients))
(define (lpf-coefficients scale)
(define s-poles (map (lambda: ([x : Complex])
(* scale x))
chebyshev-s-poles))
(define z-poles (map s-space->z-space s-poles))
(cdr (roots->coefficients z-poles)))
(: chebyshev-s-poles Poles)
(define chebyshev-s-poles
(let ()
(define num-poles 4)
(define epsilon 1.0)
(define left-half
(for/list: : (Listof Complex)
([m (in-range num-poles)])
(* i (cos (+ (* (/ 1 num-poles)
(acos (/ i epsilon)))
(/ (* pi m) num-poles))))))
left-half))
(: fir-filter ((Listof (List Nonnegative-Fixnum Real)) -> Signal -> Signal))
(define (fir-filter params)
(match params
[`((,delays ,amplitudes) ...)
(lambda: ([signal : Signal])
(define max-delay
(up-to-power-of-two (+ 1 (apply max delays))))
(: delay-buf (Vectorof Real))
(define delay-buf (make-vector max-delay 0.0))
(define next-idx 0)
(define last-t -1)
(: delays/t (Listof Nonnegative-Integer))
(define delays/t
(cond [(andmap exact-nonnegative-integer? delays)
delays]
[(error 'impossible "Make TR happy")]))
(: amplitudes/t (Listof Real))
(define amplitudes/t
(cond [(andmap inexact-real? amplitudes)
amplitudes]
[(error 'impossible "Make TR happy")]))
(lambda (t)
(unless (= t (add1 last-t))
(error
'fir-filter
"called with t=~s, expecting t=~s. Sorry about that limitation."
t
(add1 last-t)))
(let ([this-val (signal t)])
(begin
(vector-set! delay-buf next-idx this-val)
(define result
(for/fold:
([sum : Real 0.0])
([d (in-list delays/t)]
[a (in-list amplitudes/t)])
(+ sum
(* a
(vector-ref
delay-buf
(modulo (- next-idx d) max-delay))))))
(set! last-t (add1 last-t))
(set! next-idx (modulo (add1 next-idx) max-delay))
result))))]
[other (raise-type-error 'fir-filter "(listof (list number number))" 0 params)]))
(: s-space->z-space (Complex -> Complex))
(define (s-space->z-space pole) (exp pole))
(: sum-of ((Listof Complex) -> Complex))
(define (sum-of l) (foldl + 0.0+0.0i l))
(: product-of ((Listof Complex) -> Complex))
(define (product-of l) (foldl * 1.0+0.0i l))
(: all-but-n (All (T) (Natural (Listof T) -> (Listof (Listof T)))))
(define (all-but-n n l)
(cond [(= n 0) (list l)]
[(= n (length l)) (list '())]
[else (define drop-this-one
(all-but-n (- n 1) (cdr l)))
(define keep-this-one
(map (lambda: ([x : (Listof T)])
(cons (car l) x))
(all-but-n n (cdr l))))
(append drop-this-one keep-this-one)]))
(: real-part/ck (Complex -> Real))
(define (real-part/ck i)
(define angl (angle i))
(define wrapped-angle (cond [(< angl (- (/ pi 2))) (+ angl (* 2 pi))]
[else angl]))
(cond [(< (abs wrapped-angle) angle-epsilon) (magnitude i)]
[(< (abs (- pi wrapped-angle)) angle-epsilon)
(- (magnitude i))]
[else (error 'real-part/ck "angle ~s of complex number ~s is not close to zero or pi." wrapped-angle i)]))
(define angle-epsilon 1e-5)
(: up-to-power-of-two (Positive-Fixnum -> Exact-Positive-Integer))
(define (up-to-power-of-two n)
(define log-2 (log 2))
(define log-2-n (log (max n 1)))
(cond [(<= log-2-n 0) (error 'up-to-power-of-two
"impossible, make TR happy.")]
[(<= log-2 0) (error 'up-to-power-of-two
"impossible, make TR happy.")]
[else
(expt 2 (ceiling
(inexact->exact
(/ log-2-n log-2))))]))
(: my-signal Signal)
(define my-signal
(lambda (x) (sin (exact->inexact
(* 2.0 pi (/ 1.0 44100.0) 400.0 x)))))