comparison src/scheme/old/genfilter.scm @ 0:bf79fb79ee13

Initial Mercurial check in.

+;;; 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)
+	(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 z hop)
+	(on-change W (.assign K (total-matrix W f)))
+	(matexec W "load")
+)
+
+(define (ica-recon s A sink hop)
+	(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)
+
+