Skip to content

Commit

Permalink
Cleaning up Cry Baby NDK model
Browse files Browse the repository at this point in the history
  • Loading branch information
@jatinchowdhury18
jatinchowdhury18 committed Jun 16, 2023
1 parent 26af098 commit 6cc6da0
Show file tree
Hide file tree
Showing 5 changed files with 241 additions and 182 deletions.
2 changes: 2 additions & 0 deletions modules/cmake/WarningFlags.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,8 @@ elseif((CMAKE_CXX_COMPILER_ID STREQUAL "Clang") OR (CMAKE_CXX_COMPILER_ID STREQU
# For XSIMD
-Wno-cast-align -Wno-shadow -Wno-implicit-int-conversion
-Wno-zero-as-null-pointer-constant -Wno-sign-conversion
# For Eigen
-Wno-deprecated-anon-enum-enum-conversion
# Needed for ARM processor, OSX versions below 10.14
-fno-aligned-allocation
)
Expand Down
23 changes: 11 additions & 12 deletions src/processors/other/cry_baby/CryBaby.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,9 @@ void CryBaby::prepare (double sampleRate, int samplesPerBlock)
}

ndk_model = std::make_unique<CryBabyNDK>();
ndk_model->prepare ((needsOversampling ? oversampleRatio : 1.0) * sampleRate, 0.5);
ndk_model->reset ((needsOversampling ? oversampleRatio : 1.0) * sampleRate);
const auto alpha = (double) alphaSmooth.getCurrentValue();
ndk_model->update_pots ({ (1.0 - alpha) * CryBabyNDK::VR1, alpha * CryBabyNDK::VR1 });

// pre-buffering
AudioBuffer<float> buffer (2, samplesPerBlock);
Expand All @@ -79,24 +81,21 @@ void CryBaby::processAudio (AudioBuffer<float>& buffer)
if (alphaSmooth.isSmoothing())
{
const auto alphaSmoothData = alphaSmooth.getSmoothedBuffer();
for (int sample = 0; sample < block.getNumSamples();)
for (auto [ch, n, data] : chowdsp::buffer_iters::sub_blocks<32, true> (block))
{
const auto samplesToProcess = juce::jmin (32, block.getNumSamples() - sample);
if (ch == 0)
{
const auto alpha = (double) alphaSmoothData[n / smootherDivide];
ndk_model->update_pots ({ (1.0 - alpha) * CryBabyNDK::VR1, alpha * CryBabyNDK::VR1 });
}

const auto subBufferView = chowdsp::BufferView<float> { block, sample, samplesToProcess };

ndk_model->set_alpha ((double) alphaSmoothData[sample / smootherDivide]);
for (auto [ch, data] : chowdsp::buffer_iters::channels (subBufferView))
ndk_model->process_channel (data, ch);

sample += samplesToProcess;
ndk_model->process (data, ch);
}
}
else
{
ndk_model->set_alpha ((double) alphaSmooth.getCurrentValue());
for (auto [ch, channelData] : chowdsp::buffer_iters::channels (block))
ndk_model->process_channel (channelData, ch);
ndk_model->process (channelData, ch);
}
};

Expand Down
248 changes: 112 additions & 136 deletions src/processors/other/cry_baby/CryBabyNDK.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,12 @@
/*
* This file was generated on 2023-06-15 22:18:56.182984
* using the command: `generate_ndk_cpp.py cry_baby_config.json`
*/
#include "CryBabyNDK.h"

namespace
{
// START USER ENTRIES
constexpr auto R1 = 68.0e3;
constexpr auto R2 = 22.0e3;
constexpr auto R3 = 390.0;
Expand All @@ -12,7 +17,6 @@ constexpr auto R7 = 33.0e3;
constexpr auto R8 = 82.0e3;
constexpr auto R9 = 10.0e3;
constexpr auto R10 = 1.0e3;
constexpr auto VR1 = 100.0e3;
constexpr auto C1 = 10.0e-9;
constexpr auto C2 = 4.7e-6;
constexpr auto C3 = 10.0e-9;
Expand All @@ -21,138 +25,116 @@ constexpr auto C5 = 220.0e-9;
constexpr auto L1 = 500.0e-3;
constexpr auto Vcc = 9.0;
constexpr auto Vt = 26.0e-3;
constexpr auto Is = 20.3e-15;
constexpr auto BetaF = 1430.0;
constexpr auto BetaR = 4.0;
constexpr auto Is_Q1 = 20.3e-15;
constexpr auto BetaF_Q1 = 1430.0;
constexpr auto AlphaF_Q1 = (1.0 + BetaF_Q1) / BetaF_Q1;
constexpr auto BetaR_Q1 = 4.0;
constexpr auto Is_Q2 = 20.3e-15;
constexpr auto BetaF_Q2 = 1430.0;
constexpr auto AlphaF_Q2 = (1.0 + BetaF_Q2) / BetaF_Q2;
constexpr auto BetaR_Q2 = 4.0;
// END USER ENTRIES
} // namespace

void CryBabyNDK::prepare (double fs, double initial_alpha)
void CryBabyNDK::reset (T fs)
{
const auto Gr = Eigen::DiagonalMatrix<double, 10> { 1.0 / R1, 1.0 / R2, 1.0 / R3, 1.0 / R4, 1.0 / R5, 1.0 / R6, 1.0 / R7, 1.0 / R8, 1.0 / R9, 1.0 / R10 };
const auto Gx = Eigen::DiagonalMatrix<double, 6> { 2.0 * fs * C1, 2.0 * fs * C2, 2.0 * fs * C3, 2.0 * fs * C4, 2.0 * fs * C5, 1.0 / (2.0 * fs * L1) };
const auto Z = Eigen::DiagonalMatrix<double, 6> { 1.0, 1.0, 1.0, 1.0, 1.0, -1.0 };

const Eigen::Matrix<double, 10, 13> Nr {
{ +0, +0, +0, +0, +0, +0, +0, +0, -1, +0, +1, +0, +0 },
{ -1, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +1, +0 },
{ +0, +0, +0, +0, +0, +0, +0, +1, +0, +0, +0, +0, +0 },
{ +1, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, -1 },
{ -1, +0, +0, +0, +0, +0, +0, +0, +0, +1, +0, +0, +0 },
{ +0, +0, +0, +1, -1, +0, +0, +0, +0, +0, +0, +0, +0 },
{ +0, +0, +0, +0, -1, +0, +0, +0, +0, +0, +0, +0, +1 },
{ +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +1 },
{ +0, +0, +0, +0, +0, -1, +0, +0, +0, +0, +0, +0, +0 },
{ +0, +0, +0, +0, +0, +0, -1, +0, +0, +0, +0, +1, +0 },
const auto Gr = Eigen::DiagonalMatrix<T, num_resistors> { (T) 1 / R1, (T) 1 / R2, (T) 1 / R3, (T) 1 / R4, (T) 1 / R5, (T) 1 / R6, (T) 1 / R7, (T) 1 / R8, (T) 1 / R9, (T) 1 / R10 };
const auto Gx = Eigen::DiagonalMatrix<T, num_states> { (T) 2 * fs * C1, (T) 2 * fs * C2, (T) 2 * fs * C3, (T) 2 * fs * C4, (T) 2 * fs * C5, (T) 1 / ((T) 2 * fs * L1) };
const auto Z = Eigen::DiagonalMatrix<T, num_states> { 1, 1, 1, 1, 1, -1 };

// Set up component-defined matrices
const Eigen::Matrix<T, num_resistors, num_nodes> Nr {
{ +1, -1, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0 },
{ +0, +0, +0, -1, +0, +0, +0, +0, +0, +0, +0, +0, +1 },
{ +0, +0, +0, +0, +1, +0, +0, +0, +0, +0, +0, +0, +0 },
{ +0, +0, +0, +1, +0, -1, +0, +0, +0, +0, +0, +0, +0 },
{ +0, +0, +0, +1, +0, +0, +0, +0, +0, -1, +0, +0, +0 },
{ +0, +0, +1, +0, +0, +0, -1, +0, +0, +0, +0, +0, +0 },
{ +0, +0, +0, +0, +0, +1, -1, +0, +0, +0, +0, +0, +0 },
{ +0, +0, +0, +0, +0, +1, +0, +0, +0, +0, +0, +0, +0 },
{ +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +1, +0 },
{ +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, -1, +0, +1 },
};

const Eigen::Matrix<double, 2, 13> Nv {
{ +0, +0, +1, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0 },
{ +0, +1, -1, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0 },
const Eigen::Matrix<T, num_pots, num_nodes> Nv {
{ +0, +0, +0, +0, +0, +0, +0, +0, +1, +0, +0, +0, +0 },
{ +0, +0, +0, +0, +0, +0, +0, +1, -1, +0, +0, +0, +0 },
};

const Eigen::Matrix<double, 6, 13> Nx {
{ +0, +0, +0, -1, +0, +0, +0, +0, +1, +0, +0, +0, +0 },
{ +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +1 },
{ +0, +0, +0, +0, +1, -1, +0, +0, +0, +0, +0, +0, +0 },
{ -1, +1, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0 },
{ +0, +0, +1, +0, +0, +0, +0, +0, +0, -1, +0, +0, +0 },
{ +0, +0, +0, +0, -1, +0, +0, +0, +0, +0, +0, +0, +1 },
const Eigen::Matrix<T, num_states, num_nodes> Nx {
{ +0, -1, +1, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0 },
{ +0, +0, +0, +0, +0, +1, +0, +0, +0, +0, +0, +0, +0 },
{ +0, +0, +0, +0, +0, +0, +1, +0, +0, +0, +0, -1, +0 },
{ +0, +0, +0, +1, +0, +0, +0, -1, +0, +0, +0, +0, +0 },
{ +0, +0, +0, +0, +0, +0, +0, +0, +1, -1, +0, +0, +0 },
{ +0, +0, +0, +0, +0, +1, -1, +0, +0, +0, +0, +0, +0 },
};

const Eigen::Matrix<double, 2, 13> Nu {
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 },
const Eigen::Matrix<T, num_voltages, num_nodes> Nu {
{ +1, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0 },
{ +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +1 },
};

const Eigen::Matrix<double, 4, 13> Nn {
{ +1, +0, +0, -1, +0, +0, +0, +0, +0, +0, +0, +0, +0 },
{ +1, +0, +0, +0, +0, +0, +0, -1, +0, +0, +0, +0, +0 },
{ +0, +0, +0, +0, +0, +0, +1, +0, +0, -1, +0, +0, +0 },
{ +0, +0, +0, +0, +0, -1, +1, +0, +0, +0, +0, +0, +0 },
const Eigen::Matrix<T, num_nl_ports, num_nodes> Nn {
{ +0, +0, -1, +1, +0, +0, +0, +0, +0, +0, +0, +0, +0 },
{ +0, +0, +0, +1, -1, +0, +0, +0, +0, +0, +0, +0, +0 },
{ +0, +0, +0, +0, +0, +0, +0, +0, +0, -1, +1, +0, +0 },
{ +0, +0, +0, +0, +0, +0, +0, +0, +0, +0, +1, -1, +0 },
};

const Eigen::Matrix<double, 1, 13> No {
{ 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
const Eigen::Matrix<T, num_outputs, num_nodes> No {
{ +0, +0, +0, +0, +0, +0, +0, +1, +0, +0, +0, +0, +0 },
};

Eigen::Matrix<double, 6, 15> Nx_0;
Nx_0.setZero();
Nx_0.block<6, 13> (0, 0) = Nx;

Eigen::Matrix<double, 4, 15> Nn_0;
Nn_0.setZero();
Nn_0.block<4, 13> (0, 0) = Nn;

Eigen::Matrix<double, 1, 15> No_0;
No_0.setZero();
No_0.block<1, 13> (0, 0) = No;

Eigen::Matrix<double, 2, 15> Zero_I;
Zero_I.setZero();
static constexpr auto S_dim = num_nodes + num_voltages;
Eigen::Matrix<T, num_states, S_dim> Nx_0 = Eigen::Matrix<T, num_states, S_dim>::Zero();
Nx_0.block<num_states, num_nodes> (0, 0) = Nx;
Eigen::Matrix<T, num_nl_ports, S_dim> Nn_0 = Eigen::Matrix<T, num_nl_ports, S_dim>::Zero();
Nn_0.block<num_nl_ports, num_nodes> (0, 0) = Nn;
Eigen::Matrix<T, num_outputs, S_dim> No_0 = Eigen::Matrix<T, num_outputs, S_dim>::Zero();
No_0.block<num_outputs, num_nodes> (0, 0) = No;
Eigen::Matrix<T, num_pots, S_dim> Nv_0 = Eigen::Matrix<T, num_pots, S_dim>::Zero();
Nv_0.block<num_pots, num_nodes> (0, 0) = Nv;
Eigen::Matrix<T, num_voltages, S_dim> Zero_I = Eigen::Matrix<T, num_voltages, S_dim>::Zero();
Zero_I.block<2, 2> (0, 13).setIdentity();

//============================================
Eigen::Matrix<double, 15, 15> S0_mat;
S0_mat.setZero();
S0_mat.block<13, 13> (0, 0) = Nr.transpose() * Gr * Nr
+ Nx.transpose() * Gx * Nx;
S0_mat.block<13, 2> (0, 13) = Nu.transpose();
S0_mat.block<2, 13> (13, 0) = Nu;
Eigen::Matrix<double, 15, 15> S0_inv = S0_mat.inverse();

Eigen::Matrix<double, 2, 15> Nv_0;
Nv_0.setZero();
Nv_0.block<2, 13> (0, 0) = Nv;
Eigen::Matrix<T, S_dim, S_dim> S0_mat = Eigen::Matrix<T, S_dim, S_dim>::Zero();
S0_mat.block<num_nodes, num_nodes> (0, 0) = Nr.transpose() * Gr * Nr + Nx.transpose() * Gx * Nx;
S0_mat.block<num_nodes, num_voltages> (0, num_nodes) = Nu.transpose();
S0_mat.block<num_voltages, num_nodes> (num_nodes, 0) = Nu;
Eigen::Matrix<T, S_dim, S_dim> S0_inv = S0_mat.inverse();

// Pre-compute NDK and intermediate matrices
Q = Nv_0 * (S0_inv * Nv_0.transpose());
Ux = Nx_0 * (S0_inv * Nv_0.transpose());
Uo = No_0 * (S0_inv * Nv_0.transpose());
Un = Nn_0 * (S0_inv * Nv_0.transpose());
Uu = Zero_I * (S0_inv * Nv_0.transpose());

A0_mat = 2.0 * (Z * (Gx * (Nx_0 * (S0_inv * Nx_0.transpose())))) - Z.toDenseMatrix();
B0_mat = 2.0 * (Z * (Gx * (Nx_0 * (S0_inv * Zero_I.transpose()))));
C0_mat = 2.0 * (Z * (Gx * (Nx_0 * (S0_inv * Nn_0.transpose()))));
A0_mat = (T) 2 * (Z * (Gx * (Nx_0 * (S0_inv * Nx_0.transpose())))) - Z.toDenseMatrix();
B0_mat = (T) 2 * (Z * (Gx * (Nx_0 * (S0_inv * Zero_I.transpose()))));
C0_mat = (T) 2 * (Z * (Gx * (Nx_0 * (S0_inv * Nn_0.transpose()))));
D0_mat = No_0 * (S0_inv * Nx_0.transpose());
E0_mat = No_0 * (S0_inv * Zero_I.transpose());
F0_mat = No_0 * (S0_inv * Nn_0.transpose());
G0_mat = Nn_0 * (S0_inv * Nx_0.transpose());
H0_mat = Nn_0 * (S0_inv * Zero_I.transpose());
K0_mat = Nn_0 * (S0_inv * Nn_0.transpose());
two_Z_Gx = 2.0 * (Z.toDenseMatrix() * Gx.toDenseMatrix());

prev_alpha = -1.0; // force the update!
set_alpha (initial_alpha);
two_Z_Gx = (T) 2 * (Z.toDenseMatrix() * Gx.toDenseMatrix());

//============================================
// reset state
for (size_t ch = 0; ch < 2; ++ch)
// reset state vectors
for (size_t ch = 0; ch < MAX_NUM_CHANNELS; ++ch)
{
x_n[ch].setZero();
// v_n[ch].setZero();
// x_n[ch] = Eigen::Vector<double, 6> { -0.00129042068087091,
// 0.6066797237167123,
// -0.0064893267388429644,
// -0.19469253168118281,
// -0.19741507823194129,
// -1.5319085469643699E-7 };
v_n[ch] = Eigen::Vector<double, 4> { 3.9271560942528319, 4.524363916506168, 3.9262980403171812, 4.5429223634080538 };
v_n[ch] = Eigen::Vector<T, num_nl_ports> { 3.9, 4.5, 3.9, 4.5 };
}
}

void CryBabyNDK::set_alpha (double alpha) noexcept
void CryBabyNDK::update_pots (const std::array<T, num_pots>& pot_values)
{
if (alpha == prev_alpha)
return;

jassert (alpha >= 0.1 && alpha <= 0.99);
alpha = juce::jlimit (0.1, 0.99, alpha);

prev_alpha = alpha;

const auto Rv = Eigen::Matrix<double, 2, 2> { { (1.0 - alpha) * VR1, 0.0 }, { 0.0, alpha * VR1 } };

Eigen::Matrix<double, 2, 2> Rv_Q_inv = (Rv + Q).inverse();
Eigen::Matrix<T, num_pots, num_pots> Rv = Eigen::Matrix<T, num_pots, num_pots>::Zero();
Rv (0, 0) = pot_values[0];
Rv (1, 1) = pot_values[1];
Eigen::Matrix<T, num_pots, num_pots> Rv_Q_inv = (Rv + Q).inverse();

A_mat = A0_mat - (two_Z_Gx * (Ux * (Rv_Q_inv * Ux.transpose())));
B_mat = B0_mat - (two_Z_Gx * (Ux * (Rv_Q_inv * Uu.transpose())));
Expand All @@ -165,76 +147,70 @@ void CryBabyNDK::set_alpha (double alpha) noexcept
K_mat = K0_mat - (Un * (Rv_Q_inv * Un.transpose()));
}

void CryBabyNDK::process_channel (std::span<float> x, size_t ch) noexcept
void CryBabyNDK::process (std::span<float> channel_data, size_t ch) noexcept
{
Eigen::Vector<double, 2> u_n { 0.0, Vcc };

Eigen::Vector<double, 4> p_n;
Eigen::Matrix<double, 4, 4> Jac;
Jac.setZero();
Eigen::Vector<double, 4> i_n;
Eigen::Vector<double, 4> F_min;
Eigen::Matrix<double, 4, 4> A_solve;
const Eigen::Matrix<double, 4, 4> eye4 = Eigen::Matrix<double, 4, 4>::Identity();
Eigen::Vector<double, 4> delta_v;
Eigen::Vector<double, 1> y_n;

for (auto& sample : x)
Eigen::Vector<T, num_voltages> u_n { (T) 0, Vcc };
Eigen::Vector<T, num_nl_ports> p_n;
Eigen::Matrix<T, num_nl_ports, num_nl_ports> Jac = Eigen::Matrix<T, num_nl_ports, num_nl_ports>::Zero();
Eigen::Vector<T, num_nl_ports> i_n;
Eigen::Vector<T, num_nl_ports> F_min;
Eigen::Matrix<T, num_nl_ports, num_nl_ports> A_solve;
const Eigen::Matrix<T, num_nl_ports, num_nl_ports> eye = Eigen::Matrix<T, num_nl_ports, num_nl_ports>::Identity();
Eigen::Vector<T, num_nl_ports> delta_v;
Eigen::Vector<T, num_outputs> y_n;

for (auto& sample : channel_data)
{
u_n (0) = (double) sample;

u_n (0) = (T) sample;
p_n.noalias() = G_mat * x_n[ch] + H_mat * u_n;

static constexpr auto alphaF = (1 + BetaF) / BetaF;
double exp_v1_v0;
double exp_mv0;
double exp_v3_v2;
double exp_mv2;
T exp_v1_v0;
T exp_mv0;
T exp_v3_v2;
T exp_mv2;
const auto calc_currents = [&]
{
exp_v1_v0 = std::exp ((v_n[ch](1) - v_n[ch](0)) / Vt);
exp_mv0 = std::exp (-v_n[ch](0) / Vt);
i_n (0) = Is_Q1 * ((exp_v1_v0 - (T) 1) / BetaF_Q1 + (exp_mv0 - (T) 1) / BetaR_Q1);
i_n (1) = -Is_Q1 * (-(exp_mv0 - (T) 1) + AlphaF_Q1 * (exp_v1_v0 - (T) 1));

exp_v3_v2 = std::exp ((v_n[ch](3) - v_n[ch](2)) / Vt);
exp_mv2 = std::exp (-v_n[ch](2) / Vt);

i_n (0) = Is * ((exp_v1_v0 - 1.0) / BetaF + (exp_mv0 - 1.0) / BetaR);
i_n (1) = -Is * (-(exp_mv0 - 1.0) + alphaF * (exp_v1_v0 - 1.0));
i_n (2) = Is * ((exp_v3_v2 - 1.0) / BetaF + (exp_mv2 - 1.0) / BetaR);
i_n (3) = -Is * (-(exp_mv2 - 1.0) + alphaF * (exp_v3_v2 - 1.0));
i_n (2) = Is_Q2 * ((exp_v3_v2 - (T) 1) / BetaF_Q2 + (exp_mv2 - (T) 1) / BetaR_Q2);
i_n (3) = -Is_Q2 * (-(exp_mv2 - (T) 1) + AlphaF_Q2 * (exp_v3_v2 - (T) 1));
};

const auto calc_jacobian = [&]
{
Jac (0, 0) = Is / Vt * (-exp_v1_v0 / BetaF - exp_mv0 / BetaR);
Jac (0, 1) = Is / Vt * (exp_v1_v0 / BetaF);
Jac (1, 0) = Is / Vt * (-exp_mv0 + alphaF * exp_v1_v0);
Jac (1, 1) = Is / Vt * (-alphaF * exp_v1_v0);
Jac (2, 2) = Is / Vt * (-exp_v3_v2 / BetaF - exp_mv2 / BetaR);
Jac (2, 3) = Is / Vt * (exp_v3_v2 / BetaF);
Jac (3, 2) = Is / Vt * (-exp_mv2 + alphaF * exp_v3_v2);
Jac (3, 3) = Is / Vt * (-alphaF * exp_v3_v2);
Jac (0, 0) = (Is_Q1 / Vt) * (-exp_v1_v0 / BetaF_Q1 - exp_mv0 / BetaR_Q1);
Jac (0, 1) = (Is_Q1 / Vt) * (exp_v1_v0 / BetaF_Q1);
Jac (1, 0) = (Is_Q1 / Vt) * (-exp_mv0 + AlphaF_Q1 * exp_v1_v0);
Jac (1, 1) = (Is_Q1 / Vt) * (-AlphaF_Q1 * exp_v1_v0);

Jac (2, 2) = (Is_Q2 / Vt) * (-exp_v3_v2 / BetaF_Q2 - exp_mv2 / BetaR_Q2);
Jac (2, 3) = (Is_Q2 / Vt) * (exp_v3_v2 / BetaF_Q2);
Jac (3, 2) = (Is_Q2 / Vt) * (-exp_mv2 + AlphaF_Q2 * exp_v3_v2);
Jac (3, 3) = (Is_Q2 / Vt) * (-AlphaF_Q2 * exp_v3_v2);
};

double delta;
T delta;
int nIters = 0;
do
{
calc_currents();
calc_jacobian();

F_min.noalias() = p_n + K_mat * i_n - v_n[ch];
A_solve.noalias() = K_mat * Jac - eye4;
A_solve.noalias() = K_mat * Jac - eye;
delta_v.noalias() = A_solve.householderQr().solve (F_min);
v_n[ch] -= delta_v;

delta = delta_v.array().abs().sum();
} while (delta > 1.0e-2 && ++nIters < 8);

calc_currents();

y_n.noalias() = D_mat * x_n[ch] + E_mat * u_n + F_mat * i_n;
sample = (float) y_n (0);

x_n[ch] = A_mat * x_n[ch] + B_mat * u_n + C_mat * i_n;
}
}
Loading

0 comments on commit 6cc6da0

Please sign in to comment.