-
Notifications
You must be signed in to change notification settings - Fork 3
/
burg.lisp
194 lines (171 loc) · 6.32 KB
/
burg.lisp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
(defpackage :burg
(:use :cl)
(:export
:memcof
:predict
:filter
:blit
:vec
))
(in-package :burg)
;; -------------------------------------------------
;; Burg AR Analysis
(defvar *burg-library*
#+:MAC (translate-logical-pathname "PROJECTS:DYLIB;libLispBurg.dylib")
#+:WIN32 "Lisp_FPU.dll")
(fli:register-module *burg-library*)
(fli:define-foreign-function (lisp_memcof
"lisp_memcof" :source)
((src (:pointer :double))
(n :long)
(dst (:pointer :double))
(m :long))
:result-type :double
:language :ansi-c)
(defun memcof (arr m)
;; returns polynomial coeffs of order ord -> ord+1 coffs
(let* ((n (length arr))
(ans (make-array m :element-type 'double-float)))
(fli:with-dynamic-foreign-objects ()
(let ((cdst (fli:allocate-dynamic-foreign-object
:type :double
:nelems m))
(csrc (fli:allocate-dynamic-foreign-object
:type :double
:nelems n)))
(dotimes (ix n)
(setf (fli:dereference csrc :index ix)
(coerce (aref arr ix) 'double-float)))
(let ((xms (lisp_memcof csrc n cdst m)))
(dotimes (ix m)
(setf (aref ans ix) (fli:dereference cdst :index ix)))
(values xms ans))
))))
#|
(let* ((nformants 3)
(ord (* 2 nformants))
(m ord)
(err 0.250d0)
(ndata 100)
(y (vm:vectorwise ((x (vm:dramp ndata))
(noise (vm:gnoise ndata :sd err)))
(+ noise
(sin (* 2.0d0 pi x 0.1d0))
(* 0.25 (sin (* 2d0 pi x 0.25)))
(* -0.06 (sin (* 2d0 pi x 0.35)))))))
(labels ((frq (z)
(/ (phase z) 2 pi))
(bw (z)
(/ (log (abs z)) -2 2 pi)))
(multiple-value-bind (xms dcof)
(memcof y m)
(plt:plot 'coffs dcof :clear t)
(let* ((nfft 16384)
(nfft/2 (truncate nfft 2))
(arr (let ((arr (make-array nfft
:element-type 'single-float
:initial-element 0.0)))
(setf (aref arr 0) 1.0)
(loop for coff across dcof
for ix from 1
do
(setf (aref arr ix) (- (coerce coff 'single-float))))
arr))
(fcoffs (fft:fwd-magnitude-db arr))
(frqs (vops:vscale (/ nfft) (vm:framp nfft/2)))
(fcoffs (map 'vector '- fcoffs)))
(plt:plot 'dbplt frqs fcoffs :clear t
;; :xrange '(0.09 0.11)
))
(plt:plot 'data y :clear t
:yrange '(-5 5))
(multiple-value-bind (xms ym)
;; (predict y m :nnew 80 :nprev 0)
(filter y m)
(plt:plot 'data ym :color :red)
(plt:plot 'diff (map 'vector #'- ym y) :clear t)
(let* ((poly-coffs (Let ((arr (make-array (1+ ord)
:element-type 'double-float)))
(setf (aref arr 0) 1d0)
(loop for coff across dcof
for ix from 1
do
(setf (aref arr ix) (- coff)))
(nreverse arr)))
(roots (remove-if (lambda (z)
(minusp (imagpart z)))
(nrglue:zroots poly-coffs))))
(list :rms-err xms
:coffs dcof
:roots (map 'vector (lambda (z)
(complex (abs z)
(/ (phase z) 2 pi)))
roots)
:freq (map 'vector #'frq roots)
:bw (map 'vector #'bw roots))
)))))
|#
;; --------------------------------------------------------
;; Vector Routines
(defun blit (src srcoff dst dstoff nel)
(loop for ix from 0 below nel do
(setf (aref dst (+ dstoff ix))
(aref src (max 0 (+ srcoff ix))))))
(defun vec (arr offs &optional (nel (- (length arr) offs)))
(make-array nel
:displaced-to arr
:displaced-index-offset offs
:element-type (array-element-type arr)))
;; ------------------------------------------------------
;; Burg AR Prediction
(defun predict (data m
&key
nnew nprev
(mean-remove t))
(let* ((nel (length data))
(mn (if mean-remove (vm:mean data) 0.0d0))
(dat (if mean-remove
(vm:vectorwise ((d data))
(- d mn))
data)))
(multiple-value-bind (xms dcof) (memcof dat m)
(let* ((dx (nreverse dcof))
(ans (make-array (+ m nprev nnew))))
(blit dat (- nel (+ nprev m)) ans 0 (+ nprev m))
(loop for kx from (+ nprev m) below (+ nprev m nnew) do
(setf (aref ans kx)
(vm:inner-prod dx (vec ans (- kx m) m))
))
(loop for kx from (+ nprev m -1) downto m do
(setf (aref ans kx)
(vm:inner-prod dx (vec ans (- kx m) m))
))
(values xms
(vm:vectorwise ((a (subseq ans m)))
(+ a mn)))
))))
(defun filter (data m
&key
(mean-remove t))
(let* ((nel (length data))
(mn (if mean-remove (vm:mean data) 0.0d0))
(dat (if mean-remove
(vm:vectorwise ((d data))
(- d mn))
data)))
(multiple-value-bind (xms dcof) (memcof dat m)
(let* ((dx (nreverse dcof))
(ans (make-array nel)))
(loop for pos from 0 below nel do
(setf (aref ans pos)
(+ mn ;; (aref data pos)
(loop for ix from (- pos m) below pos
for coff across dx
sum
(if (minusp ix)
0
(* coff (aref data ix)))
))
))
(values xms ans)
))))