Improve making roots legal by increasing lower stat weights in fit

See #16.

src/psd_model_fitting.jl

@@ -37,19 +37,30 @@ function fit_psd(PSD::AbstractVector, pulsemodel::AbstractVector, p, q=-1)
     w = pulseFT2 ./ PSD.^3
     # Don't let the DC bin have ZERO weight, else model likes to go negative, particularly
     # If there's lots of "action" (poles) near ω=0, or cos(ω)=1.
-    w[1] = w[2]*.01
-
-    aaa_hybrid = aaawt(z, PSD, w, p)
-    vfit = vectorfit(z, PSD, w, aaa_hybrid.poles, q)
-
-    vfit = make_poles_legal(vfit, z, PSD, w)
-    vfit, cosroots = make_roots_legal(vfit)
-
-    # zpoles = exp.(acosh.(vfit.λ))
-    # zroots = exp.(acosh.(complex(cosroots)))
-    # var = mean(PSD)
-    model = ARMAModel(vfit)
-    vfit, model
+    w[1] = w[2]*.2
+
+    # Signal models tend to put very little weight at ω=π (like, 1e-9x the max weight).
+    # This can lead to terrible fits up there, as well as illegal roots. Prevent both by checking
+    # 1. ...for illegal roots
+    # 2. ...for minimum function value < 0.1x the minimum target function value.
+    # If either are found, boost the minimum weight and repeat.
+    minpsd = minimum(PSD)
+    while true
+        aaa_hybrid = aaawt(z, PSD, w, p)
+        vfit = vectorfit(z, PSD, w, aaa_hybrid.poles, q)
+
+        vfit = make_poles_legal(vfit, z, PSD, w)
+        minmodel = minimum(real(vfit(z)))
+        if minmodel > 0.2*minpsd
+            cosroots = RCPRoots(roots(vfit))
+            illegal_roots = illegal_RPs(cosroots)
+            if length(illegal_roots) == 0
+                model = ARMAModel(vfit)
+                return vfit, model
+            end
+        end
+        w .+= 2minimum(w)  # Triple min weight each iteration
+    end
 end
 
 """
@@ -80,25 +91,25 @@ constant. Other strategies for that case are TBD.
 """
 function make_roots_legal(vfit::PartialFracRational)
     while true
-        ma_roots = RCPRoots(roots(vfit))
-        illegal_roots = illegal_RPs(ma_roots)
+        cosroots = RCPRoots(roots(vfit))
+        illegal_roots = illegal_RPs(cosroots)
         if length(illegal_roots) == 0
-            return vfit, ma_roots
+            return vfit, cosroots
         end
-        # println("Have to fix $(length(illegal_roots)) illegal roots")
+        # println("Have to fix $(length(illegal_roots)) illegal roots: $illegal_roots")
 
         if vfit.m ≥ vfit.n
             # Strategy, when m≥n: add a constant to vfit until roots are legal.
             r = illegal_roots
             midpts = 0.5*(r[1:end-1] .+ r[2:end])
             testpts = vcat(midpts, [-1,0,1])
             f = real(vfit(testpts))
-            b = vfit.b
+            b = copy(vfit.b)
             b[1] -= 1.01*minimum(f)
             vfit = PartialFracRational(vfit.λ, vfit.a, b)
         else
             # Strategy, when m<n: I don't have one!
-            @show ma_roots
+            @show cosroots
             @show vfit
             throw(ErrorException("No strategy for making roots legal for m<n."))
         end

src/vector_fitting.jl

@@ -31,7 +31,7 @@
 # n(z) function will contain all the information, the d(z) denominator will be unity, and the nonlinear
 # coefficients b will all be zero. In that limit, n(z) and R(z) will contain all the remaining unknowns, and
 # they can be found in a linear problem. That is, fitting F to data samples s becomes equivalent to
-# fitting F(Z)*d(z)=n(z)+d(z)*R(z) to s*d(z) in the limit that d → 1. We try to approach that limit
+# fitting F(z)*d(z)=n(z)+d(z)*R(z) to s*d(z) in the limit that d → 1. We try to approach that limit
 # iteratively.
 #
 # Based on Gustavsen, B., & Semlyen, A. (1999). "Rational approximation of frequency domain responses by
@@ -127,10 +127,10 @@ function vectorfit(z::AbstractVector, f::AbstractVector, wt::AbstractVector, λ0
     end
     WM = W*[C Φ]
     optparam = WM\Wf
-    model = real([C Φ]*optparam)
-
     ρ = optparam[1:n]
     c = real(optparam[n+1:end])
+
+    # model = real([C Φ]*optparam)
     # Wtres = norm(W*model.-Wf)
     # @show Wtres, ρ, c