;nyquist plug-in
;version 4
;type generate
;name "VFD..."
;preview linear
;action "Generating VFD..."
;author "Duncan Yang"
;release "3.2.0"
;copyright "Released under terms of the GNU General Public License version 2 or later."

;control carrier "Carrier frequency" float-text "Hz" 1000.0 1.0 1000000.0
;control f "Modulator frequency" float-text "Hz" 60.0 0.1 200000.0
;control dur "Duration" time "" 1 0 3600
;control mp "Mod depth" float "%" 90 -100 100
;control dt "Dead time" float "μs" 0.0 0.0 50.0
;control h3inj "Third harmonic injection" choice "off,on" 0
;control levels "Output levels" choice "2-level,3-level,5-level" 0
;control bias "Width offset" float "%" 0 0 100
;control amp "Amplitude" float "%" 100 0 100
;control smooth "Smoothing" choice "Enabled,Disabled" 0

;; --- Convert controls to internal units ---
(setq mp (* 0.01 mp))
(setq bias (* 0.01 bias))
(setq dt (* 1.0e-6 dt))

(defun build-modulator (inject-h3 hz duration phase)
    (if (= inject-h3 1)
        (scale (/ 2.0 (sqrt 3.0))
            (sum (osc (hz-to-step hz) duration *table* phase)
                (scale (/ 1.0 6.0)
                    (osc (hz-to-step (* 3.0 hz)) duration *table* (* 3.0 phase)))))
        (osc (hz-to-step hz) duration *table* phase)))

(defun apply-dt (binary t-dead gen-dur)
;; Idealized dead-time: during each switching transition the leg dwells at the
;; midpoint (average of the two rails) for t-dead -- both switches off, with no
;; commitment to a current direction. output(t) = (binary(t) + binary(t-dt)) / 2.
;; This is a linear, symmetric operation: it trims |output| equally at every
;; edge, producing no current-dependent crossover distortion (and therefore
;; needs no compensation). Input gen-dur long, output (gen-dur - t-dead) long.
    (if (> t-dead 0)
        (let* ((head (extract-abs 0 (- gen-dur t-dead) (cue binary)))
                (delayed (seq (s-rest t-dead) head))
                (avg (scale 0.5 (sum binary delayed))))
            (extract-abs t-dead gen-dur (cue avg)))
        binary))

(defun tri-wave (frq dur)
;; Sharp +/-1 triangle carrier with NO wavetable oscillator: osc band-limits
;; every table in this build, rounding any non-sine waveform toward a sine,
;; and a sinusoidal carrier warps the sine reference into a triangle via the
;; arcsin mapping. Instead, take the sign of a cosine (a sharp square via clip)
;; and integrate it to get genuine linear ramps, scale to +/-1, then high-pass
;; to remove any integrator DC drift. Only osc(sine), clip, integrate, scale,
;; and hp are used -- all of which behave correctly here.
    (hp (scale (* 4.0 frq)
            (integrate (clip (scale 1.0e6 (osc (hz-to-step frq) dur *table* 90.0)) 1.0)))
        (/ frq 100.0)))

(defun pwm-2level (ref tri t-dead gen-dur)
    (apply-dt (clip (scale 1.0e6 (diff ref tri)) 1.0) t-dead gen-dur))

(defun pwm-3level (ref tri t-dead gen-dur)
;; Single H-bridge unipolar PWM. Two legs compared to +ref and -ref, each given
;; the idealized dwell-at-0 dead-time, then bridged.
    (let* ((ref-neg (scale -1.0 ref))
            (leg-a (apply-dt (clip (scale 1.0e6 (diff ref tri)) 1.0) t-dead gen-dur))
            (leg-b (apply-dt (clip (scale 1.0e6 (diff ref-neg tri)) 1.0) t-dead gen-dur)))
        (scale 0.5 (diff leg-a leg-b))))

(defun pwm-5level (ref tri frq t-dead gen-dur)
;; Cascaded H-bridge PS-PWM. Two bridges in series, each producing 3-level,
;; carriers shifted by T_c/4. Each leg (4 total) gets the idealized dwell-at-0
;; dead-time before bridging.
    (let* ((tc-quarter (/ 1.0 (* 4.0 frq)))
            (tri2 (extract-abs tc-quarter (+ tc-quarter gen-dur)
                (cue (tri-wave frq (+ gen-dur tc-quarter)))))
            (ref-neg (scale -1.0 ref))
            (a1 (apply-dt (clip (scale 1.0e6 (diff ref tri)) 1.0) t-dead gen-dur))
            (b1 (apply-dt (clip (scale 1.0e6 (diff ref-neg tri)) 1.0) t-dead gen-dur))
            (a2 (apply-dt (clip (scale 1.0e6 (diff ref tri2)) 1.0) t-dead gen-dur))
            (b2 (apply-dt (clip (scale 1.0e6 (diff ref-neg tri2)) 1.0) t-dead gen-dur))
            (bridge1 (scale 0.5 (diff a1 b1)))
            (bridge2 (scale 0.5 (diff a2 b2))))
        (scale 0.5 (sum bridge1 bridge2))))

(defun pwm (frq duration width modamp modfrq inject-h3 levels-mode smooth-on t-dead)
    (let* ((target-sr *sound-srate*)
            (oversample 16)
            (high-sr (* target-sr oversample))
            (anti-alias (* 0.4999999 target-sr))
            (wl (if smooth-on (/ (float (+ 1 (truncate (* 0.001 frq)))) frq) 0.0))
            (mod-phase (- (* 360.0 modfrq (+ wl t-dead))))
            (gen-dur (+ duration wl t-dead))
            (wave (progv '(*sound-srate*) (list high-sr)
                (let* ((m (build-modulator inject-h3 modfrq gen-dur mod-phase))
                        (ref (sum width (scale modamp m)))
                        (tri (tri-wave frq gen-dur))
                        (pwm (cond
                            ((= levels-mode 2) (pwm-5level ref tri frq t-dead gen-dur))
                            ((= levels-mode 1) (pwm-3level ref tri t-dead gen-dur))
                            (t (pwm-2level ref tri t-dead gen-dur)))))
                    (force-srate target-sr (lowpass8 pwm anti-alias))))))
        (when smooth-on
            (setq wave (lowpass2 wave (/ *sound-srate* 3.8) 0.55))
            (setq wave (extract-abs wl (+ duration wl) (cue wave))))
        wave))

(cond
    ((> carrier (/ *sound-srate* 2)) ;; Enforce: Carrier must be less than 1/2 the sample rate
        (print "Error:\n\nCarrier frequency will alias at current sample rate."))
    ((> (* 5 f) carrier) ;; Enforce: Carrier must also be at least 5x the modulator frequency
        (print "Error:\n\nCarrier frequency must be at least 5 times the modulator frequency."))
    ((> dt (/ 0.1 carrier)) ;; Make sure dead time is not excessive
        (print "Error:\n\nDead time cannot take up more than 1/10 of each period."))
    ((> mp (- 0.99 (* dt carrier)))
        (print "Error:\n\nModulation index too high for dead time."))
    ((<= amp 0) ;; Check for sub-zero amplitude
        (print "Error:\n\nAmplitude is non-positive, nothing to do."))
    ((<= dur 0) ;; Check for zero duration
        (print "Error:\n\nDuration is non-positive, nothing to do."))
    ;; Else generate PWM
    (T (scale (* 0.01 amp)
        (pwm carrier dur bias mp f h3inj levels (= smooth 0) dt))))