samer@0: ;;; Preamble
samer@0: (load "audio.scm")
samer@0: (load "lineout.scm")
samer@0: (load "functions.scm")
samer@0: 
samer@0: (import "samer.mds.*")
samer@0: 
samer@0: ;;; transfer filter fn to vector of values
samer@0: (define (set-filter F fn)
samer@0: 	(define f (.mat F))
samer@0: 	(for-each
samer@0: 		(lambda (i) (.set f 0 i (fn i)))
samer@0: 		(range 0 N))
samer@0: 	(.changed F)
samer@0: )
samer@0: 
samer@0: ;; convert vector of values into a filter fn
samer@0: (define (filter-fn F)
samer@0: 	(let ((mat (.mat F)))
samer@0: 		(lambda (i) (.get mat 0 i))))
samer@0: 
samer@0: ;; returns fn f(i) that returns 1 if d'th coordinate of P(i) <op> th
samer@0: ;; eg (coor-test P 2 > 0.5) is function which returns P(i,2) > 0.5
samer@0: (define (coor-test P d op th)
samer@0: 	(lambda (i) (if (op (.get P i d) th) 1.0 0.0)))
samer@0: 
samer@0: 
samer@0: ;;; set up tasks for ICA domain filtering using f as a function
samer@0: ;;; which returns (f i) as the scaling factor for the ith ICA component
samer@0: ;;; source is an AudioSource.
samer@0: (define (direct-recon source N hop f sink)
samer@0: 	(define x (linein source N hop))
samer@0: 	(define W (Matrix. "W" N N))		; ICA matrix
samer@0: 	(define K (Matrix. "K" N N))			; total filtering matrix
samer@0: 	(define z (VVector. "z" N))				; reconstruction
samer@0: 	(addtask (Ops.times z K x) (Ops.update z))
samer@0: 	(overlap-and-add sink z hop)
samer@0: 	(on-change W (.assign K (total-matrix W f)))
samer@0: 	(matexec W "load")
samer@0: )
samer@0: 
samer@0: (define (ica-recon s A sink hop)
samer@0: 	(define N (.size s))
samer@0: 	(define F (VVector. "F" N))			; the source domain filter
samer@0: 	(define z (VVector. "z" N))
samer@0: 	(addtasks (Ops.timesEquals s F) (Ops.times z A s))
samer@0: 	(overlap-and-add sink z hop))
samer@0: 
samer@0: 
samer@0: (define (recon source N hop sink)
samer@0: 	(define x (norm (linein source N hop)))
samer@0: 	(define ica (mk-ica x laplace-spec))
samer@0: 	(define s (.output ica))
samer@0: 	(define F (VVector. "F" N))			; the source domain filter
samer@0: 	(define W (.getWeightMatrix ica))
samer@0: 	(define A (.getBasisMatrix ica))
samer@0: 	(define z (VVector. "z" N))
samer@0: 	(define order (Matrix. "order" 5 512))
samer@0: 
samer@0: 	(addtasks (Ops.timesEquals s F) (Ops.times z A s))
samer@0: 	(overlap-and-add sink z hop)
samer@0: 
samer@0: 	(on-change W (exec ica "basis"))
samer@0: 	(matexec W "load")
samer@0: )
samer@0: 
samer@0: ;; returns total transformation matrix  formed by filtering
samer@0: ;; in the ICA domain. f is a function which returns the scaling
samer@0: ;; factor for the i'th component, ie (f i) is a scalar.
samer@0: (define (total-matrix W f)
samer@0: 	(define D (Jama.Matrix. N N 0.0))
samer@0: 	(for-each
samer@0: 		(lambda (i) (.set D i i (f i)))
samer@0:       (range 0 N))
samer@0: 	(.times (.inverse W) (.times D W))
samer@0: )
samer@0: 
samer@0: ; eg (total-matrix W (coor-test P 1 > 0.0))
samer@0: ; filters out all components whos MDS y-coor <= 0
samer@0: 
samer@0: ;; transfer values from G to F reordered by order
samer@0: (define-method (rank-filter F order G)
samer@0: 	(rank-filter F order G 0 (.size G)))
samer@0: 
samer@0: (define-method (rank-filter F order G start end)
samer@0: 	(Mathx.set F (Zero.))
samer@0: 	(let ((f (.mat F)) (g (.mat G)) (I (.mat order)))
samer@0: 		(for-each
samer@0: 			(lambda (i) (.set f 0  (.get I 0 i)  (.get g 0 i)))
samer@0: 			(range start end)))
samer@0: 	(.changed F)
samer@0: )
samer@0: 
samer@0: ;;; filter is a certain spherically symmetric kernel centred on
samer@0: ;;; a particular component. Geometry directly from correlation
samer@0: ;;; matrix (not from MDS configuation)
samer@0: (define (proximity-filter R F fn)
samer@0: 	(ProximityFilter. R F (VFunction. "proximity kernel" fn)) f)
samer@0: 
samer@0: (define (geometric-filter P f)
samer@0: 	(define n (.getRowDimension P))
samer@0: 	(define E (.getColumnDimension P))
samer@0: 	(define off (VDouble. "offset" 0.0))
samer@0: 	(define norm	(VVector. "normal" E))
samer@0: 	(define hp	(LogisticHyperplane. (.array norm) (.value$ off)))
samer@0: 	(define gf		(GeometricFilter. P f hp))
samer@0: 	(define obs (observer (.setOffset hp (.value$ off)) (.run gf) (.changed f)))
samer@0: 	(.addObserver P obs)
samer@0: 	(.addObserver off obs)
samer@0: 	(.addObserver norm obs)
samer@0: 	f)
samer@0: 
samer@0: