# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

import sys
from __future__ import division

# +
import numpy as np
from phasor.utilities.ipynb.displays import *
from phasor.utilities.ipynb.sympy import *

import declarative

from declarative.bunch import (
    DeepBunch
)

import phasor.math.dispatched as dmath 
#import phasor.math.dispatch_sympy

# +
import phasor.utilities.version as version
print(version.foundations_version())

from phasor.utilities.np import logspaced

from phasor import optics
from phasor import base 
from phasor import signals 
from phasor import system
from phasor import readouts 

from phasor.signals.pade_fit.pade_fit import ratsvd

import scipy.signal

# +
F_nyquist = 16384 / 2
F_AC = F_nyquist * 2 * np.arange(0, 1001) / 301
    
F_AC = logspaced(1, 7000, 500)
ZPKz = (
    [
        .9,
        -.3,
        .1+.95j,
        .1-.95j,
    ], 
    [
        #.9,
        .99,
        #.97,
        -0.03+.95j,
        -.03-.95j,
        #-.2,
    ], 
    10
)


b,a = scipy.signal.zpk2tf(*ZPKz)


Fb = mplfigB(Nrows=2)
w, h = scipy.signal.freqz_zpk(*ZPKz , worN = F_AC / F_nyquist * np.pi)
F_exact = h
Fb.ax0.loglog(F_AC, abs(h), label = 'discrete')
Fb.ax1.semilogx(F_AC, np.angle(h), label = 'discrete')

w, h = scipy.signal.freqz(b, a , worN = F_AC / F_nyquist * np.pi)
#F_exact = h
Fb.ax0.loglog(F_AC, abs(h), label = 'discrete')
Fb.ax1.semilogx(F_AC, np.angle(h), label = 'discrete')

Fb.ax0.legend()
# -

b

R = np.random.randn(len(F_AC)) + 1j*np.random.randn(len(F_AC))
F_noise = F_exact * (1 + R / 10)

fitz = ratsvd(F_AC, F_noise, W = 10, 
              order_a = 10, order_b = 10, 
              F_nyquist= F_nyquist,
              max_size = 100,
              ms_method = 'sum',
              #aorder_min = 300,
             )
#fits

# +
Fb = mplfigB(Nrows=2)
w, h = scipy.signal.freqz(b, a, worN = F_AC / F_nyquist * np.pi)
Fb.ax0.loglog(F_AC, abs(h), label = 'discrete')
Fb.ax1.semilogx(F_AC, np.angle(h), label = 'discrete')

h = F_exact
Fb.ax0.loglog(F_AC, abs(F_exact), label = 'Fex')
Fb.ax1.semilogx(F_AC, np.angle(F_exact), label = 'Fex')

h = F_exact
Fb.ax0.loglog(F_AC, abs(F_noise), label = 'Fnoise')
Fb.ax1.semilogx(F_AC, np.angle(F_noise), label = 'Fnoise')


for ftup in fitz[:4]:
    res, b_fit, a_fit, order_b, order_a, F_rescale = ftup
    w, h = scipy.signal.freqz(b_fit, a_fit, worN = F_AC / F_rescale * np.pi)
    Fb.ax0.loglog(
        F_AC, 
        abs(h), 
        label = 'fit {order_b} ({res})'.format(
            order_a = order_a, 
            order_b = order_b,
            res = res,
        ),
        #color = 'green',
    )
    Fb.ax1.semilogx(F_AC, np.angle(h), label = 'fit')

Fb.ax0.legend(ncol = 3, fontsize = 6)
# -

fits = ratsvd(F_AC, F_exact, W = 10, 
              order_a = 10, order_b = 10, 
              F_nyquist=None,
              max_size = 100,
              ms_method = 'sum',
              #aorder_min = 300,
             )
fits

# +
Fb = mplfigB(Nrows=2)
w, h = scipy.signal.freqz(b, a, worN = F_AC / F_nyquist * np.pi)
Fb.ax0.loglog(F_AC, abs(h), label = 'discrete')
Fb.ax1.semilogx(F_AC, np.angle(h), label = 'discrete')

h = F_exact
Fb.ax0.loglog(F_AC, abs(F_exact), label = 'Fex')
Fb.ax1.semilogx(F_AC, np.angle(F_exact), label = 'Fex')

h = F_exact
Fb.ax0.loglog(F_AC, abs(F_noise), label = 'Fnoise')
Fb.ax1.semilogx(F_AC, np.angle(F_noise), label = 'Fnoise')



for ftup in fits[:4]:
    res, b_fit, a_fit, order_b, order_a, F_scale = ftup
    w, h = scipy.signal.freqs(b_fit, a_fit, worN = F_AC / F_rescale * 2 * np.pi)
    Fb.ax0.loglog(
        F_AC, 
        abs(h), 
        label = 'fit {order_b} ({res})'.format(
            order_a = order_a, 
            order_b = order_b,
            res = res,
        ),
        #color = 'green',
    )
    Fb.ax1.semilogx(F_AC, np.angle(h), label = 'fit')

Fb.ax0.legend(ncol = 3, fontsize = 6)
# -





StarcoderdataJupyterScriptsDedupFiltered