diff options
author | 2020-05-03 17:23:50 +0200 | |
---|---|---|
committer | 2020-05-03 17:23:50 +0200 | |
commit | b2c40086e6b003121e61719da93944f4e360e33d (patch) | |
tree | 47e7c257323a89c2e9a529f2c1817fd52d723a50 /src | |
parent | b5be886e01cd04c410783997dcf8469e1a2733a2 (diff) |
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.
Diffstat (limited to 'src')
-rw-r--r-- | src/lib/netlist/devices/net_lib.cpp | 2 | ||||
-rw-r--r-- | src/lib/netlist/devices/nld_system.cpp | 9 | ||||
-rw-r--r-- | src/lib/netlist/devices/nld_system.h | 17 | ||||
-rw-r--r-- | src/lib/netlist/devices/nlid_system.h | 60 | ||||
-rw-r--r-- | src/lib/netlist/examples/noise.cpp | 36 | ||||
-rw-r--r-- | src/lib/netlist/nl_base.h | 11 | ||||
-rw-r--r-- | src/lib/netlist/plib/pmath.h | 9 | ||||
-rw-r--r-- | src/lib/netlist/plib/prandom.h | 136 |
8 files changed, 255 insertions, 25 deletions
diff --git a/src/lib/netlist/devices/net_lib.cpp b/src/lib/netlist/devices/net_lib.cpp index 7e23ecfa2f6..2ed1733b041 100644 --- a/src/lib/netlist/devices/net_lib.cpp +++ b/src/lib/netlist/devices/net_lib.cpp @@ -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) diff --git a/src/lib/netlist/devices/nld_system.cpp b/src/lib/netlist/devices/nld_system.cpp index 1888c1b9fb0..b5689d637f5 100644 --- a/src/lib/netlist/devices/nld_system.cpp +++ b/src/lib/netlist/devices/nld_system.cpp @@ -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", "") diff --git a/src/lib/netlist/devices/nld_system.h b/src/lib/netlist/devices/nld_system.h index 879361f3ee1..9da5b32890d 100644 --- a/src/lib/netlist/devices/nld_system.h +++ b/src/lib/netlist/devices/nld_system.h @@ -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) \ diff --git a/src/lib/netlist/devices/nlid_system.h b/src/lib/netlist/devices/nlid_system.h index b318ad05491..2e55d181e9f 100644 --- a/src/lib/netlist/devices/nlid_system.h +++ b/src/lib/netlist/devices/nlid_system.h @@ -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 diff --git a/src/lib/netlist/examples/noise.cpp b/src/lib/netlist/examples/noise.cpp new file mode 100644 index 00000000000..870d526e477 --- /dev/null +++ b/src/lib/netlist/examples/noise.cpp @@ -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] diff --git a/src/lib/netlist/nl_base.h b/src/lib/netlist/nl_base.h index 08edcf9307e..a504dc800ba 100644 --- a/src/lib/netlist/nl_base.h +++ b/src/lib/netlist/nl_base.h @@ -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) { diff --git a/src/lib/netlist/plib/pmath.h b/src/lib/netlist/plib/pmath.h index cb6ae6fc8cd..5c5f8704716 100644 --- a/src/lib/netlist/plib/pmath.h +++ b/src/lib/netlist/plib/pmath.h @@ -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 diff --git a/src/lib/netlist/plib/prandom.h b/src/lib/netlist/plib/prandom.h index 5e26ef4997d..5aebb9ccca4 100644 --- a/src/lib/netlist/plib/prandom.h +++ b/src/lib/netlist/plib/prandom.h @@ -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::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<T, N> mt; - std::size_t index; + std::array<FT, 256> m_buf; + std::size_t m_p; + FT m_stddev; }; using mt19937_64 = mersenne_twister_t< |