Refactor clipfix.ny to fix multiple problems

Basically a rewrite based on the original algorithm.

plug-ins/clipfix.ny

@@ -1,141 +1,100 @@
 ;nyquist plug-in
-;version 1
+;version 4
 ;type process
 ;preview enabled
-;categories "http://audacityteam.org/namespace#NoiseRemoval"
 ;name "Clip Fix..."
 ;action "Reconstructing clips..."
-;author "Benjamin Schwartz"
+;author "Benjamin Schwartz and Steve Daulton"
 ;copyright "Licensing confirmed under terms of the GNU General Public License version 2"
 
-;; clipfix.ny by Benjamin Schwartz.
+;; Algorithm by Benjamin Schwartz
-;; Licensing confirmed under terms of the GNU General Public License version 2:
+;; Clip Fix is a simple, stupid (but not blind) digital-clipping-corrector
-;; http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
+;; The algorithm is fairly simple:
-;; with kind agreement of Benjamin Schwartz, December 2011.
+;; 1. Find all clipped regions
-;; GUI updated by Steve Daulton July 2012
+;; 2. Get the slope immediately on either side of the region
-;;
+;; 3. Do a cubic spline interpolation.
-;; For information about writing and modifying Nyquist plug-ins:
+;; 4. Go to next region
-;; http://wiki.audacityteam.org/wiki/Nyquist_Plug-ins_Reference
 
-;control thresh "Threshold of Clipping (%)" real "" 95 0 100
+;control threshold "Threshold of Clipping (%)" float "" 95 0 100
+;control gain "Reduce amplitude to allow for restored peaks (dB)" float "" -9 -30 0
 
-(setf largenumber 100000000) ;;Largest number of samples that can be imported
+(setf threshold (/ threshold 100))
-(setf blocksize 100000)
+(setf gain (db-to-linear gain))
-
+(setf buffersize 100000)
-;;Clip Fix is a simple, stupid (but not blind) digital-clipping-corrector
+(setf slopelength 4)  ; number of samples used to calculate the exit / re-entry slope
-;;The algorithm is fairly simple:
-;;1. Find all clipped regions
-;;2. Get the slope immediately on either side of the region
-;;3. Do a cubic spline interpolation.
-;;4. Go to next region
-
-;;Coded from start (didn't know lisp (well, scheme, but not not lisp and certainly not
-;;some XLISP 2.0 derivative)) to finish
-;;(fully working, more or less) in one afternoon (and some evening).
-;;Written by Benjamin Schwartz, MIT class of 2006, on May 25, 2004.
-;;Explanatory text added by Gale Andrews, May 2008.
-
-(defun declip (sin) ;;Central function
-(let* ((threshold  (* (peak sin largenumber) thresh 0.01))
-(s2 (snd-copy sin))
-(samplerate (snd-srate s2))
-(s2length (snd-length s2 largenumber)))
-
-(seqrep (i (1+ (/ s2length blocksize)))
-  (let ((l (min blocksize (- s2length (* i blocksize)))))
-     ;;(print (list i t0 l samplerate))
-     (snd-from-array 0 samplerate 
-	(workhorse 
-		;;(let () (print (list s2 (type-of s2) l (type-of l)))
-			(snd-fetch-array s2 l l)
-		;;)
-			threshold))))
-
-;;(setf r (snd-fetch-array (snd-copy s) (snd-length s largenumber) 1)) ;;Create a sound array
-;;(snd-from-array (snd-t0 s) (snd-srate s) (workhorse r threshold))
-))
-
-(defun workhorse (r threshold)
-
-(setf n (length r)) ;; Record its length
-
-(setf exithigh ()) ;;Times when the wavefrom left the allowed region
-(setf returnhigh ())  ;;Times when it returned to the allowed region
-
-(setf drange 4)
-
-(let ((i drange) (max (- n drange))) ;;Leave room at ends for derivative processing
-  (while (< i max)
-  (if (>= (aref r i) threshold) 
-    (if (< (aref r (- i 1)) threshold)
-      (setq exithigh (cons (- i 1) exithigh))) ;;We just crossed the threshold up
-    (if (>= (aref r (- i 1)) threshold)
-      (setq returnhigh (cons i returnhigh)))) ;;We just crossed the threshold down
-  (setq i (1+ i))))
-
-(setq exithigh (reverse exithigh)) ;;List comes out backwards
-(setq returnhigh (reverse returnhigh))
-
-(if (>= (aref r (1- drange)) threshold) ;;If the audio begins in a clipped region, ignore 
-  (setq returnhigh (cdr returnhigh))) ;the extra return from threshold
-
-(setf exitlow ())  ;; Same as above, but for the bottom threshold
-(setf returnlow ())
-
-(setf threshlow (* -1 threshold)) ;;Assumes your digital range is zero-centered
 
 
-(let ((i drange) (max (- n drange)))
+(defun declip (sig thresh peak)
-  (while (< i max)
+  (let* ((threshold (* thresh peak))
-  (if (<= (aref r i) threshlow)
+         (ln (truncate len))
-    (if (> (aref r (- i 1)) threshlow)
+         (finalbufsize (rem ln buffersize)))
-      (setq exitlow (cons (- i 1) exitlow)))
+    ;; Calculate the number of buffers we can process.
-    (if (<= (aref r (- i 1)) threshlow)
+    ;; if final buffer is not large enough for de-clipping we
-      (setq returnlow (cons i returnlow))))
+    ;; will just add it on the end as is.
-  (setq i (1+ i))))
+    (if (>= finalbufsize slopelength)
+        (setf buffercount (1+ (/ ln buffersize))) 
+        (setf buffercount (/ ln buffersize)))
+    ;;; Make output sequence from processed buffers
+    (setf out
+      (seqrep (i buffercount)
+        (let* ((step (min buffersize (- ln (* i buffersize))))
+               (buffer (snd-fetch-array sig step step))
+               (processed (process buffer threshold step)))
+          (cue (mult gain
+                    (snd-from-array 0 *sound-srate* processed))))))
+    ;;; If there's unprocessed audio remaining, add it to the end
+    (if (and (> finalbufsize 0)(< finalbufsize slopelength))
+        (seq out (cue (getfinalblock sig finalbufsize gain)))
+        out)))
 
-(setq exitlow (reverse exitlow))
-(setq returnlow (reverse returnlow))
 
-(if (<= (aref r (1- drange)) threshlow)
+(defun getfinalblock (sig step gain)
-  (setq returnlow (cdr returnlow)))
+  (let ((block (snd-fetch-array sig step step)))
+    (mult gain (snd-from-array 0 *sound-srate* block))))
 
-(while (and exithigh returnhigh) ;;If there are more clipped regions
-    (let* ((t1 (car exithigh))  ;;exit time
-          (t2 (car returnhigh)) ;;return time
-          (d1 (max 0 (/ (- (aref r t1) (aref r (- t1 (1- drange)))) (1- drange)))) ;;slope at exit
-          (d2 (min 0 (/ (- (aref r (+ t2 (1- drange))) (aref r t2)) (1- drange)))) ;;slope at return
-          (m (/ (+ d2 d1) (* (- t2 t1) (- t2 t1)))) ;;interpolation is by (t-t1)(t-t2)(mx+b)
-          (b (- (/ d2 (- t2 t1)) (* m t2))) ;;These values of m and b make the cubic seamless
-          (j (1+ t1))) ;; j is the index
 
-     (while (< j t2)
+(defun process (buffer threshold bufferlength)
-         (setf (aref r j) (+ (aref r t1) (* (- j t1) (- j t2) (+ (* m j) b)))) 
+  ;;; Find threshold crossings
-	 (setf (aref r j) (+ (* (- t2 j) (/ (aref r t1) (- t2 t1))) (* (- j t1) (/ (aref r t2) (- t2 t1)))  (* (- j t1) (- j t2) (+ (* m j) b))))
+  (setf exit-list ())         ; list of times when waveform exceeds threshold
-         (setq j (1+ j)))) 
+  (setf return-list ())       ; list of times when waveform returns below threshold
-     (setq exithigh (cdr exithigh))
+  ;; Limitation of algorithm: the first and last 'slopelength' at ends of buffer are ignored
-     (setq returnhigh (cdr returnhigh)))
+  ;; so that we have enough samples beyond the threshold crossing to calculate the slope.
+  (let ((last-sample (- bufferlength slopelength)))
+    (do ((i slopelength (1+ i)))
+        ((>= i last-sample))
+      (if (>= (abs (aref buffer i)) threshold)
+          (when (< (abs (aref buffer (- i 1))) threshold)   ; we just crossed threshold
+            (push (- i 1) exit-list))
+          (when (>= (abs (aref buffer (- i 1))) threshold)  ; we just got back in range
+            (push i return-list)))))
+  ;; Reverse lists back into chronological order.
+  ;; This is faster than appending values in chronological order.
+  (setf exit-list (reverse exit-list))
+  (setf return-list (reverse return-list))
+  ;; If the audio begins in a clipped region, discard the first return
+  (when (>= (abs (aref buffer (1- slopelength))) threshold)
+    (setq return-list (cdr return-list)))
+  ;; Interpolate between each pair of exit / entry points
+  (let ((slopelen (1- slopelength)))
+    (mapc (lambda (t0 t1)
+            (interpolate buffer t0 t1 slopelen))
+          exit-list return-list))
+  buffer)
 
-(while (and exitlow returnlow) ;;Same for bottom
-    (let* ((t1 (car exitlow))
-          (t2 (car returnlow))
-          (d1 (min 0 (/ (- (aref r t1) (aref r (- t1 (1- drange)))) (1- drange)))) ;;slope at exit
-          (d2 (max 0 (/ (- (aref r (+ t2 (1- drange))) (aref r t2)) (1- drange)))) ;;slope at return
-          (m (/ (+ d2 d1) (* (- t2 t1) (- t2 t1))))
-          (b (- (/ d2 (- t2 t1)) (* m t2)))
-	  (a (/ (+ (aref r t1) (aref r t2)) 2))
-          (j (1+ t1)))
-     (while (< j t2)
-         (setf (aref r j) (+ (* (- t2 j) (/ (aref r t1) (- t2 t1))) (* (- j t1) (/ (aref r t2) (- t2 t1)))  (* (- j t1) (- j t2) (+ (* m j) b))))
-         (setq j (1+ j))))
-     (setq exitlow (cdr exitlow))
-     (setq returnlow (cdr returnlow)))
 
-r)
+(defun interpolate (buffer t0 t1 dur)
+  "Cubic spline interpolation"
+  (let* ((d0 (/ (- (aref buffer t0) (aref buffer (- t0 dur))) dur)) ; slope at start
+         (d1 (/ (- (aref buffer (+ t1 dur)) (aref buffer t1)) dur)) ; slope at end
+         (m (/ (+ d1 d0) (* (- t1 t0) (- t1 t0))))
+         (b (- (/ d1 (- t1 t0)) (* m t1))))
+    (do ((j (1+ t0) (1+ j)))
+        ((= j t1))
+      (setf (aref buffer j)
+        (+ (* (- t1 j) (/ (aref buffer t0) (- t1 t0)))
+           (* (- j t0) (/ (aref buffer t1) (- t1 t0)))
+           (* (- j t0) (- j t1) (+ (* m j) b)))))))
 
-(if (arrayp s)
-  (dotimes (j (length s))
-	(setf (aref s j) (declip (aref s j))))
-  (setq s (declip s)))
 
-s
+;; (get '*selection* 'peak) introduced in Audacity 2.1.3
+(multichan-expand #'declip *track* threshold (get '*selection* 'peak))

src/effects/nyquist/Nyquist.cpp

@@ -998,7 +998,6 @@ bool NyquistEffect::ProcessOne()
          cmd += wxString::Format(wxT("(putprop '*SELECTION* (vector %s) 'PEAK)\n"), peakString) :
          cmd += wxString::Format(wxT("(putprop '*SELECTION* %s 'PEAK)\n"), peakString);
 
-      // TODO: Documen, PEAK-LEVEL is deprecated as of 2.1.3.
       // TODO: Document, PEAK-LEVEL is nil if NaN or INF.
       if (!std::isinf(maxPeakLevel) && !std::isnan(maxPeakLevel) && (maxPeakLevel < FLT_MAX)) {
          cmd += wxString::Format(wxT("(putprop '*SELECTION* (float %s) 'PEAK-LEVEL)\n"),