netlist: Add two noise sources. [Couriersud]

The two sources act as voltage sources, though noise may also be
injected as conductivy or current noise.

SYS_NOISE_MT_U: Mersenne Twister uniform noise
SYS_NOISE_MT_N: Mersenne Twister normal noise

nld_sys_noise is templated:

	using NETLIB_NAME(sys_noise_mt_u) =
		NETLIB_NAME(sys_noise)<plib::mt19937_64,
plib::uniform_distribution_t>;

Thus the approach is scalable. The implementation is state save aware,
and thus reproducible results are guaranteed.

An example use case is provided as well, see examples/noise.cpp.
This commit is contained in:
couriersud 2020-05-03 17:23:50 +02:00
parent b5be886e01
commit b2c40086e6
9 changed files with 256 additions and 25 deletions

View File

@ -71,6 +71,7 @@ project "netlist"
MAME_DIR .. "src/lib/netlist/plib/ppmf.h",
MAME_DIR .. "src/lib/netlist/plib/ppreprocessor.cpp",
MAME_DIR .. "src/lib/netlist/plib/ppreprocessor.h",
MAME_DIR .. "src/lib/netlist/plib/prandom.h",
MAME_DIR .. "src/lib/netlist/plib/pstate.h",
MAME_DIR .. "src/lib/netlist/plib/pstonum.h",
MAME_DIR .. "src/lib/netlist/plib/pstring.cpp",

View File

@ -64,6 +64,8 @@ namespace devices
LIB_ENTRY(sys_dsw1)
LIB_ENTRY(sys_dsw2)
LIB_ENTRY(sys_compd)
LIB_ENTRY(sys_noise_mt_u)
LIB_ENTRY(sys_noise_mt_n)
LIB_ENTRY(switch1)
LIB_ENTRY(switch2)
LIB_ENTRY(nicRSFF)

View File

@ -71,6 +71,15 @@ namespace devices
NETLIB_DEVICE_IMPL(sys_dsw1, "SYS_DSW", "+I,+1,+2")
NETLIB_DEVICE_IMPL(sys_dsw2, "SYS_DSW2", "")
NETLIB_DEVICE_IMPL(sys_compd, "SYS_COMPD", "")
using NETLIB_NAME(sys_noise_mt_u) =
NETLIB_NAME(sys_noise)<plib::mt19937_64, plib::uniform_distribution_t>;
NETLIB_DEVICE_IMPL(sys_noise_mt_u, "SYS_NOISE_MT_U", "SIGMA")
using NETLIB_NAME(sys_noise_mt_n) =
NETLIB_NAME(sys_noise)<plib::mt19937_64, plib::normal_distribution_t>;
NETLIB_DEVICE_IMPL(sys_noise_mt_n, "SYS_NOISE_MT_N", "SIGMA")
NETLIB_DEVICE_IMPL(mainclock, "MAINCLOCK", "FREQ")
NETLIB_DEVICE_IMPL(gnd, "GNDA", "")
NETLIB_DEVICE_IMPL(netlistparams, "PARAMETER", "")

View File

@ -63,17 +63,20 @@
NET_C(cOUT, name.Q)
// FIXME ... remove parameters
#define SYS_DSW(name, cIN, cP1, cP2) \
NET_REGISTER_DEV(SYS_DSW, name) \
NET_C(cIN, name.I) \
NET_C(cP1, name.1) \
NET_C(cP2, name.2)
#define SYS_DSW(name, pI, p1, p2) \
NET_REGISTER_DEVEXT(SYS_DSW, name, pI, p1, p2)
#define SYS_DSW2(name) \
NET_REGISTER_DEV(SYS_DSW2, name)
NET_REGISTER_DEVEXT(SYS_DSW2, name)
#define SYS_COMPD(name) \
NET_REGISTER_DEV(SYS_COMPD, name)
NET_REGISTER_DEVEXT(SYS_COMPD, name)
#define SYS_NOISE_MT_U(name, pSIGMA) \
NET_REGISTER_DEVEXT(SYS_NOISE_MT_U, name, pSIGMA)
#define SYS_NOISE_MT_N(name, pSIGMA) \
NET_REGISTER_DEVEXT(SYS_NOISE_MT_N, name, pSIGMA)
/* Default device to hold netlist parameters */
#define PARAMETERS(name) \

View File

@ -11,6 +11,9 @@
#include "netlist/nl_base.h"
#include "netlist/nl_setup.h"
#include "netlist/plib/putil.h"
#include "netlist/plib/prandom.h"
#include <random>
namespace netlist
{
@ -520,7 +523,8 @@ namespace devices
NETLIB_UPDATEI()
{
const netlist_sig_t state = m_I();
if (1 || (state != m_last_state))
// FIXME: update should only be called if m_I changes
if (true || (state != m_last_state))
{
m_last_state = state;
//printf("Here %d\n", state);
@ -616,6 +620,60 @@ namespace devices
state_var<netlist_sig_t> m_last_state;
};
// -----------------------------------------------------------------------------
// nld_sys_noise - noise source
//
// An externally clocked noise source.
// -----------------------------------------------------------------------------
template <typename E, template<class> class D>
NETLIB_OBJECT(sys_noise)
{
public:
using engine = E;
using distribution = D<nl_fptype>;
NETLIB_CONSTRUCTOR(sys_noise)
, m_T(*this, "m_T")
, m_I(*this, "I")
, m_RI(*this, "RI", nlconst::magic(0.1))
, m_sigma(*this, "SIGMA", nlconst::zero())
, m_dis(m_sigma())
{
m_mt.save_state(*this, "m_mt");
m_dis.save_state(*this, "m_dis");
register_subalias("1", m_T.P());
register_subalias("2", m_T.N());
}
protected:
NETLIB_UPDATEI()
{
nl_fptype val = m_dis(m_mt);
m_T.change_state([this, val]()
{
m_T.set_G_V_I(plib::reciprocal(m_RI()), val, nlconst::zero());;
});
}
NETLIB_RESETI()
{
m_T.set_G_V_I(plib::reciprocal(m_RI()), nlconst::zero(), nlconst::zero());
}
private:
analog::NETLIB_SUB(twoterm) m_T;
logic_input_t m_I;
param_fp_t m_RI;
param_fp_t m_sigma;
engine m_mt;
distribution m_dis;
//std::normal_distribution<nl_fptype> m_dis;
};
} // namespace devices
} // namespace netlist

View File

@ -0,0 +1,36 @@
// license:CC0
// copyright-holders:Couriersud
/*
* noise.cpp
*
*/
//! [noise_example]
#include "netlist/devices/net_lib.h"
// ./nltool -t 1 -l X.3 -l X.4 -n oscillator src/lib/netlist/examples/noise.cpp
// X.3 : Square out
// X.4 : Triangle out
NETLIST_START(noise)
SOLVER(Solver, 48000)
CLOCK(nclk, 2000)
SYS_NOISE_MT_N(noise, 0.5)
RES(R1,1000)
RES(R2,1000)
ANALOG_INPUT(VP, 12.0)
NET_C(nclk.Q, noise.I)
NET_C(VP, R1.1, nclk.VCC)
NET_C(noise.1, R1.2)
NET_C(noise.2, R2.1)
NET_C(GND, R2.2, nclk.GND)
NETLIST_END()
//! [noise_example]

View File

@ -525,6 +525,11 @@ namespace netlist
netlist_t & exec() noexcept { return m_netlist; }
const netlist_t & exec() const noexcept { return m_netlist; }
// State saving support - allows this to be passed
// to generic save_state members
template<typename C>
void save_item(C &state, const pstring &membername, const pstring &itemname);
private:
netlist_t & m_netlist;
@ -1782,6 +1787,12 @@ namespace netlist
return m_netlist.nlstate();
}
template<typename C>
inline void detail::netlist_object_t::save_item(C &state, const pstring &membername, const pstring &itemname)
{
this->state().save(*this, state, name(), membername + "." + itemname);
}
template<class C, typename... Args>
void device_t::create_and_register_subdevice(const pstring &name, unique_pool_ptr<C> &dev, Args&&... args)
{

View File

@ -39,11 +39,16 @@ namespace plib
static inline constexpr T three() noexcept { return static_cast<T>(3); } // NOLINT
static inline constexpr T four() noexcept { return static_cast<T>(4); } // NOLINT
static inline constexpr T hundred()noexcept { return static_cast<T>(100); } // NOLINT
static inline constexpr T sqrt2() noexcept { return static_cast<T>(1.414213562373095048801688724209L); } // NOLINT
static inline constexpr T pi() noexcept { return static_cast<T>(3.14159265358979323846264338327950L); } // NOLINT
static inline constexpr T one_thirds() noexcept { return fraction(one(), three()); }
static inline constexpr T two_thirds() noexcept { return fraction(two(), three()); }
static inline constexpr T ln2() noexcept { return static_cast<T>(0.6931471805599453094172321214581766L); } // NOLINT
static inline constexpr T sqrt2() noexcept { return static_cast<T>(1.4142135623730950488016887242096982L); } // NOLINT
static inline constexpr T sqrt3() noexcept { return static_cast<T>(1.7320508075688772935274463415058723L); } // NOLINT
static inline constexpr T sqrt3_2() noexcept { return static_cast<T>(0.8660254037844386467637231707529362L); } // NOLINT
static inline constexpr T pi() noexcept { return static_cast<T>(3.1415926535897932384626433832795029L); } // NOLINT
/// \brief Electric constant of vacuum
///
static inline constexpr T eps_0() noexcept { return static_cast<T>(8.854187817e-12); } // NOLINT

View File

@ -42,8 +42,9 @@ namespace plib
class mersenne_twister_t
{
public:
mersenne_twister_t()
: index(N+1)
: m_p(N)
{
seed(5489);
}
@ -51,22 +52,29 @@ namespace plib
static constexpr T min() noexcept { return static_cast<T>(0); }
static constexpr T max() noexcept { return ~T(0) >> (sizeof(T)*8 - w); }
template <typename ST, typename STR>
void save_state(ST &st, const STR &name)
{
st.save_item(m_p, name, "index");
st.save_item(m_mt, name, "mt");
}
void seed(T val) noexcept
{
const T lowest_w(~T(0) >> (sizeof(T)*8 - w));
index = N;
mt[0] = val;
m_p = N;
m_mt[0] = val;
for (std::size_t i=1; i< N; i++)
mt[i] = (f * (mt[i-1] ^ (mt[i-1] >> (w-2))) + i) & lowest_w;
m_mt[i] = (f * (m_mt[i-1] ^ (m_mt[i-1] >> (w-2))) + i) & lowest_w;
}
T operator()() noexcept
{
const T lowest_w(~T(0) >> (sizeof(T)*8 - w));
if (index >= N)
if (m_p >= N)
twist();
T y = mt[index++];
T y = m_mt[m_p++];
y = y ^ ((y >> u) & d);
y = y ^ ((y << s) & b);
y = y ^ ((y << t) & c);
@ -77,9 +85,9 @@ namespace plib
void discard(std::size_t v) noexcept
{
if (v > N - index)
if (v > N - m_p)
{
v -= N - index;
v -= N - m_p;
twist();
}
while (v > N)
@ -87,7 +95,7 @@ namespace plib
v -= N;
twist();
}
index += v;
m_p += v;
}
private:
@ -99,15 +107,113 @@ namespace plib
for (std::size_t i=0; i<N; i++)
{
const T x((mt[i] & upper_mask) + (mt[(i+1) % N] & lower_mask));
const T x((m_mt[i] & upper_mask) + (m_mt[(i+1) % N] & lower_mask));
const T xA((x >> 1) ^ ((x & 1) ? a : 0));
mt[i] = mt[(i + m) % N] ^ xA;
m_mt[i] = m_mt[(i + m) % N] ^ xA;
}
index = 0;
}
m_p = 0;
}
std::array<T, N> mt;
std::size_t index;
std::size_t m_p;
std::array<T, N> m_mt;
};
template <typename FT, typename T>
FT normalize_uniform(T &p, FT m = constants<FT>::one(), FT b = constants<FT>::zero())
{
const auto mmin(static_cast<FT>(p.min()));
const auto mmax(static_cast<FT>(p.max()));
// -> 0 to a
return (p() - mmin) / (mmax - mmin) * m - b;
}
template<typename FT>
class uniform_distribution_t
{
public:
uniform_distribution_t(FT dev)
: m_stddev(dev) { }
template <typename P>
FT operator()(P &p) noexcept
{
// get -1 to 1
return normalize_uniform(p, constants<FT>::two(), constants<FT>::one())
* constants<FT>::sqrt3() * m_stddev;
}
template <typename ST, typename STR>
void save_state(ST &st, const STR &name)
{
/* no state to save */
}
private:
FT m_stddev;
};
template<typename FT>
class normal_distribution_t
{
public:
normal_distribution_t(FT dev)
: m_p(m_buf.size()), m_stddev(dev) { }
// Donald Knuth, Algorithm P (Polar method)
template <typename P>
FT operator()(P &p) noexcept
{
if (m_p >= m_buf.size())
fill(p);
return m_buf[m_p++];
}
template <typename ST, typename STR>
void save_state(ST &st, const STR &name)
{
st.save_item(m_p, name, "m_p");
st.save_item(m_buf, name, "m_buf");
}
private:
template <typename P>
void fill(P &p) noexcept
{
for (std::size_t i = 0; i < m_buf.size(); i += 2)
{
FT s;
FT v1;
FT v2;
do
{
v1 = normalize_uniform(p, constants<FT>::two(), constants<FT>::one()); // [-1..1[
v2 = normalize_uniform(p, constants<FT>::two(), constants<FT>::one()); // [-1..1[
s = v1 * v1 + v2 * v2;
} while (s >= constants<FT>::one());
if (s == constants<FT>::zero())
{
m_buf[i] = s;
m_buf[i+1] = s;
}
else
{
// last value with error for log(s)/s
// double: 1.000000e-305
// float: 9.999999e-37
// FIXME: with 128 bit randoms log(s)/w will fail 1/(2^128) ~ 2.9e-39
const auto m(m_stddev * plib::sqrt(-constants<FT>::two() * plib::log(s)/s));
m_buf[i] = /*mean+*/ m * v1;
m_buf[i+1] = /*mean+*/ m * v2;
}
}
m_p = 0;
}
std::array<FT, 256> m_buf;
std::size_t m_p;
FT m_stddev;
};
using mt19937_64 = mersenne_twister_t<