view src/scheme/genfilter.scm @ 1:5df24c91468d

Oh my what a mess.

;;; Preamble
(load "audio.scm")
(load "lineout.scm")
(load "functions.scm")

(import "samer.mds.*")

;;; transfer filter fn to vector of values
(define (set-filter F fn)
	(define f (.mat F))
	(for-each
		(lambda (i) (.set f 0 i (fn i)))
		(range 0 N))
	(.changed F)
)

;; convert vector of values into a filter fn
(define (filter-fn F)
	(let ((mat (.mat F)))
		(lambda (i) (.get mat 0 i))))

;; returns fn f(i) that returns 1 if d'th coordinate of P(i) <op> th
;; eg (coor-test P 2 > 0.5) is function which returns P(i,2) > 0.5
(define (coor-test P d op th)
	(lambda (i) (if (op (.get P i d) th) 1.0 0.0)))


;;; set up tasks for ICA domain filtering using f as a function
;;; which returns (f i) as the scaling factor for the ith ICA component
;;; source is an AudioSource.
(define (direct-recon source N hop f sink)
	(define x (linein source N hop))
	(define W (Matrix. "W" N N))		; ICA matrix
	(define K (Matrix. "K" N N))			; total filtering matrix
	(define z (VVector. "z" N))				; reconstruction
	(addtask (Ops.times z K x) (Ops.update z))
	(overlap-and-add sink z hop)
	(on-change W (.assign K (total-matrix W f)))
	(matexec W "load")
)

(define (ica-recon s A sink hop)
	(define N (.size s))
	(define F (VVector. "F" N))			; the source domain filter
	(define z (VVector. "z" N))
	(addtasks (Ops.timesEquals s F) (Ops.times z A s))
	(overlap-and-add sink z hop))


(define (recon source N hop sink)
	(define x (norm (linein source N hop)))
	(define ica (mk-ica x laplace-spec))
	(define s (.output ica))
	(define F (VVector. "F" N))			; the source domain filter
	(define W (.getWeightMatrix ica))
	(define A (.getBasisMatrix ica))
	(define z (VVector. "z" N))
	(define order (Matrix. "order" 5 512))

	(addtasks (Ops.timesEquals s F) (Ops.times z A s))
	(overlap-and-add sink z hop)

	(on-change W (exec ica "basis"))
	(matexec W "load")
)

;; returns total transformation matrix  formed by filtering
;; in the ICA domain. f is a function which returns the scaling
;; factor for the i'th component, ie (f i) is a scalar.
(define (total-matrix W f)
	(define D (Jama.Matrix. N N 0.0))
	(for-each
		(lambda (i) (.set D i i (f i)))
      (range 0 N))
	(.times (.inverse W) (.times D W))
)

; eg (total-matrix W (coor-test P 1 > 0.0))
; filters out all components whos MDS y-coor <= 0

;; transfer values from G to F reordered by order
(define-method (rank-filter F order G)
	(rank-filter F order G 0 (.size G)))

(define-method (rank-filter F order G start end)
	(Mathx.set F (Zero.))
	(let ((f (.mat F)) (g (.mat G)) (I (.mat order)))
		(for-each
			(lambda (i) (.set f 0  (.get I 0 i)  (.get g 0 i)))
			(range start end)))
	(.changed F)
)

;;; filter is a certain spherically symmetric kernel centred on
;;; a particular component. Geometry directly from correlation
;;; matrix (not from MDS configuation)
(define (proximity-filter R F fn)
	(ProximityFilter. R F (VFunction. "proximity kernel" fn)) f)

(define (geometric-filter P f)
	(define n (.getRowDimension P))
	(define E (.getColumnDimension P))
	(define off (VDouble. "offset" 0.0))
	(define norm	(VVector. "normal" E))
	(define hp	(LogisticHyperplane. (.array norm) (.value$ off)))
	(define gf		(GeometricFilter. P f hp))
	(define obs (observer (.setOffset hp (.value$ off)) (.run gf) (.changed f)))
	(.addObserver P obs)
	(.addObserver off obs)
	(.addObserver norm obs)
	f)