summaryrefslogtreecommitdiffstatshomepage
path: root/src
diff options
context:
space:
mode:
author couriersud <couriersud@gmx.org>2020-05-03 17:23:50 +0200
committer couriersud <couriersud@gmx.org>2020-05-03 17:23:50 +0200
commitb2c40086e6b003121e61719da93944f4e360e33d (patch)
tree47e7c257323a89c2e9a529f2c1817fd52d723a50 /src
parentb5be886e01cd04c410783997dcf8469e1a2733a2 (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.cpp2
-rw-r--r--src/lib/netlist/devices/nld_system.cpp9
-rw-r--r--src/lib/netlist/devices/nld_system.h17
-rw-r--r--src/lib/netlist/devices/nlid_system.h60
-rw-r--r--src/lib/netlist/examples/noise.cpp36
-rw-r--r--src/lib/netlist/nl_base.h11
-rw-r--r--src/lib/netlist/plib/pmath.h9
-rw-r--r--src/lib/netlist/plib/prandom.h136
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<