;nyquist plug-in
;version 1
;type process
;categories "http://audacityteam.org/namespace#NoiseRemoval"
;name "ESP Clip Fix..."
;action "Reconstructing clips..."
;info "Designed and implemented by Benjamin Schwartz.\nModified by Steve Daulton\n\nClip Fix attempts to reconstruct clipped regions by interpolating the\nlost signal. Before use, reduce amplification by 10 dB to give room for\nthe reconstruction. 'Threshold' is how close to the maximum sample\nmagnitude any sample must be to be considered clipped. If processing\nis slow, select only a few seconds of clipped audio at a time."    
;control thresh "Threshold of Clipping [%]" real "" 95 0 100
(setf largenumber 100000000) ;;Largest number of samples that can be imported
(setf blocksize 100000)

;;Clip Fix is a simple, stupid (but not blind) digital-clipping-corrector
;;The algorithm is fairly simple:
;;1. Correct any inverted peaks
;;2. Find all clipped regions
;;3. Get the slope immediately on either side of the region
;;4. Do a cubic spline interpolation.
;;5. Go to next region

;;Written by Benjamin Schwartz, MIT class of 2006, on May 25, 2004.
;;Explanatory text added by Gale Andrews, May 2008.
;;Correction of inverted peaks added by Steve Daulton Dec 2008.
;;If you modify this code, please retain the original credit to Benjamin Schwartz.


(setq slength (snd-length s largenumber))
(setq s1 s)
(setq dif 0)
(setq limit 0.5)

(setq array (snd-fetch-array s1 slength slength)) ;; create array and read in the samples
(setq current (aref array 0))
(do ((count 0 (+ count 1)))
((>= count slength) array );test expression / return expression
(setq next (aref array count))
(setq dif (diff current next))
(if (< dif 0)(setq dif (diff 0 dif)))
(if (> dif limit) (setf (aref array count) current))
(if (> dif limit) (setq next current)(setq current next))
)
(setq s (snd-from-array 0 (snd-srate s) array))

(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)))
  (while (< i max)
  (if (<= (aref r i) threshlow)
    (if (> (aref r (- i 1)) threshlow)
      (setq exitlow (cons (- i 1) exitlow)))
    (if (<= (aref r (- i 1)) threshlow)
      (setq returnlow (cons i returnlow))))
  (setq i (1+ i))))

(setq exitlow (reverse exitlow))
(setq returnlow (reverse returnlow))

(if (<= (aref r (1- drange)) threshlow)
  (setq returnlow (cdr returnlow)))

(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)
         (setf (aref r j) (+ (aref r t1) (* (- j t1) (- j t2) (+ (* m j) b)))) 
	 (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 exithigh (cdr exithigh))
     (setq returnhigh (cdr returnhigh)))

(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)

(if (arrayp s)
  (dotimes (j (length s))
	(setf (aref s j) (declip (aref s j))))
  (setq s (declip s)))

s