;nyquist plug-in
;version 1
;type analyze
;categories "http://lv2plug.in/ns/lv2core#AnalyserPlugin"
;name "Find Phonemes..."
;action "Finding sound..."
;control control-minimum-ms "Label regions not shorter than" real "ms" 20 0 1000
;control control-percentage "... in which the" real "-th" 75 0 100
;control control-frequency "... percentile of the power spectrum lies above a multiple of the frequency" real "Hz" 2000.0 0.1 10000
;control control-upto-frequency "... up to" real "Hz" 10000 0.1 22050
;control control-window-length "Window Length (samples)" int "" 256 8 32768
;control control-skip-length "Skip Length (cycles)" int "" 32 1 32768
;control control-window-choice "Window Type" choice "Hann,Rectangular" 0
;;unused
;;control control-amplitude "Amplitude (dB)" real "" -24 -60 0
;;; reverse-append list onto list
(defun revappend (xs ys)
(if (null xs)
ys
(revappend (rest xs) (cons (first xs) ys))))
;;; reverse-append array of lists into one list
(defun list-array-to-list (arr)
(let ((len (length arr)) result)
(dotimes (ii len result) (setq result (revappend (aref arr ii) result)))))
;;; last arg is a function taking a frame
(defun scan-frames (snd length-samples skip-samples window frame-visitor)
(let ((my-s (snd-copy snd)))
(do ((n 0 (1+ n))
(frame (snd-fft my-s length-samples skip-samples window)
(snd-fft my-s length-samples skip-samples window)))
((not (arrayp frame)))
(funcall frame-visitor frame))))
(defun eliminate-short-label (label-list minimum-seconds)
(if (or (null label-list)
(>= (- (cadar label-list) (caar label-list)) minimum-seconds))
label-list
(rest label-list)))
(defun push-label (label-list minimum-seconds
start end
length-samples skip-samples srate
name)
;; find the times in seconds that delimit the start of the start-th frame
;; and the end of the end-th frame
;; do not make touching or overlapping labels
(let* ((sample-period (/ 1.0 srate))
(scaled-start (* sample-period start skip-samples))
(scaled-end (* sample-period (+ length-samples (* end skip-samples)))))
(cond
((or (null label-list) (> scaled-start (cadar label-list)))
;; new label
;; decide now whether to eliminate the PREVIOUS label,
;; but not the new one
(cons (list scaled-start scaled-end name)
(eliminate-short-label label-list minimum-seconds)))
(T
;; overlap or touch the previous label, so extend it
(cons (list (caar label-list) scaled-end
name)
(rest label-list))))))
;;; should-label is a function returning the number (0 based) of a bucket
;;; or -1
;;; names is an array determining the number of buckets
(defun label-selected-frames (snd
minimum-seconds length-samples skip-samples
window should-label
names)
(let* ((srate (snd-srate snd))
(nbuckets (length names))
(results (make-array nbuckets)) ; array of lists
(starts-of-labels (make-array nbuckets))
(nn 0)
(make-labels
#'(lambda (frame)
(let ((bucket (funcall should-label frame)))
(dotimes (ii (1+ bucket))
(if (not (aref starts-of-labels ii))
;; transition in
(setf (aref starts-of-labels ii) nn)))
(do ((ii (1+ bucket) (1+ ii)))
((= ii nbuckets))
(if (numberp (aref starts-of-labels ii))
;; transition out
(setf (aref results ii)
(push-label (aref results ii)
minimum-seconds
(aref starts-of-labels ii) (1- nn)
length-samples skip-samples srate
(aref names ii))
(aref starts-of-labels ii) nil))))
(setq nn (1+ nn)))))
(scan-frames snd length-samples skip-samples window make-labels)
(dotimes (ii 0 nbuckets)
(let ((start-of-label (aref starts-of-labels ii)))
(setf (aref results ii)
(if (not start-of-label)
(eliminate-short-label (aref results ii) minimum-seconds)
(push-label (aref results ii)
minimum-seconds start-of-label (1- nn)
length-samples skip-samples srate
(aref names ii))))))
results))
;;; takes an array, returns an array.
;;; We will examine odd-numbered members of the output and ignore the evens.
;;; Array member 1 will be proportional to energy in DC band,
;;; 3 to energy in DC plus first sine and cosine, 5 to that plus second, etc.
(defun cumul-power-spectrum (frame)
;; make dc comparable to the other components
(setf (aref frame 0) (* 0.5 (aref frame 0)))
(let* ((snd (snd-from-array 0.0 1.0 frame)) ; the array as a sound
(sndsq (mult snd snd)) ; the squares
;; following finds partial sums, always gives 0 in first
;; place, and ignores the last component of sndsq which is
;; the Nyquist frequency. Would divide it by 2 first as with DC
;; and pad with a zero if I cared.
(integral (integrate sndsq))
(len (length frame)))
;; be nice and leave no net side effect.
(setf (aref frame 0) (* 2 (aref frame 0)))
;; convert sound back to an array.
(snd-fetch-array integral len len)))
(defun scan-power-spectra (snd
minimum-seconds fraction
;;abovebelow
freq-hz
upto-freq-hz
;;amp-db
length-samples skip-samples
window-choice)
(let* ((srate (snd-srate snd))
;;stuff from a different experiment
;; (length-samples (truncate (* length-cycles period-samples)))
;; (skip-samples (truncate (* skip-cycles period-samples)))
(window-length-seconds (/ length-samples srate))
(window-freq-hz (/ window-length-seconds))
(window
(case window-choice
(0 (mult 0.5
(sum 1 (osc (hz-to-step window-freq-hz)
window-length-seconds *sine-table* -90))))
(1 nil)))
;; ignore this stuff from a different experiment
;; ;; take the threshold in dB, convert to linear, and multiply
;; ;; by length-samples / 2 to make it comparable
;; ;; to the coefficients computed
;; ;; by fft which have the same factor in them (apart from dc and
;; ;; Nyquist frequencies for which it is length-samples).
;; (amp-linear (db-to-linear amp-db))
;; (amp-linear-scaled (* amp-linear length-samples 0.5))
;; (ampsq (* amp-linear-scaled amp-linear-scaled))
;; (my-fun
;; #'(lambda (frame)
;; (let* ((coeff-odd (aref frame index-odd))
;; (coeff-even (aref frame index-even))
;; (coeffsq (+ (* coeff-odd coeff-odd)
;; (* coeff-even coeff-even))))
;; (>= coeffsq ampsq))))
(nbuckets (truncate (/ upto-freq-hz freq-hz)))
(my-fun2
#' (lambda (frame)
(let* ((len (length frame))
(sums (cumul-power-spectrum frame))
(nsums (length sums))
(sum (aref sums (1- nsums)))
(threshold (* fraction sum))
;; scan sums for first ODD entry at or above threshold,
;; then return the EVEN index into frame for the
;; coefficient pair for that frequency
(index (do ((n 0 (+ 2 n)))
((or (>= (1+ n) nsums)
(>= (aref sums (1+ n)) threshold)) n)))
;; find the frequency corresponding to index
(threshold-frequency (* window-freq-hz (/ index 2)))
;; what multiple of the chosen frequency?
(ratio (/ threshold-frequency freq-hz)))
(min (1- nbuckets) (1- (truncate ratio))))))
(names (make-array nbuckets)))
;; initialize names
(dotimes (ii nbuckets) (setf (aref names ii)
(format nil "~A" (* (1+ ii) freq-hz))))
(label-selected-frames snd
minimum-seconds length-samples skip-samples window
my-fun2
names)))
(or
(list-array-to-list
(scan-power-spectra s
(/ control-minimum-ms 1000.0) (/ control-percentage 100.0)
;;control-abovebelow
control-frequency
control-upto-frequency
;; control-amplitude
control-window-length control-skip-length
control-window-choice))
"Found nothing")