6997134815

Committed 26 Nov 2023 05:57PM UTC coverage: 74.663% (-10.1%) from 84.742%

Build # 6997134815

Build Type

push

github

Committed by

aewallin

Commit Message

coverage workflow

Run Details

1273 of 1705 relevant lines covered (74.66%)

2.17 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

97.78

/allantools/realtime.py

"""
    Real-time statistics
    ====================

    **Author:** Anders Wallin (anders.e.e.wallin "at" gmail.com)

    Statistics are computed 'on-the-fly' from a stream of
    phase/frequency samples
    Initial version 2019 July, Anders E. E. Wallin.

    This file is part of allantools, see https://github.com/aewallin/allantools

    Version history
    ---------------

    **2019.07**
    - Initial release

    License
    -------
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.
    You should have received a copy of the GNU Lesser General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""

import numpy


class dev_realtime(object):
    """ Base-class for real-time statistics """
    def __init__(self, afs=[1], tau0=1.0, auto_afs=False, pts_per_decade=4):
        self.x = []                         # phase time-series
        self.afs = afs                      # averaging factor, tau = af*tau0
        self.auto_afs = auto_afs
        # logscpace will fail at >=6 (?), need to remove duplicates?
        self.af_taus = numpy.logspace(0, 1, pts_per_decade+1)[:-1]
        # logspace index, to keep track of afs in auto-af mode
        self.af_idx = 0
        self.af_decade = 0    # logspace decade
        if auto_afs:
            self.afs = numpy.array([1])
        self.dev = numpy.zeros(len(afs))    # resulting xDEV
        self.tau0 = tau0                     # time-interval between points

    def update_af(self):
        """ used in auto-AF mode,
            - check if we can add another AF
            - if yes, add it.
        """
        next_idx = self.af_idx+1
        next_decade = self.af_decade
        if next_idx == len(self.af_taus):
            next_idx = 0
            next_decade = self.af_decade + 1

        # next possible AF:
        next_af = int(numpy.round(
            pow(10.0, next_decade) * self.af_taus[next_idx]))
        if len(self.x) >= (2*next_af+1):  # can compute next AF
            self.afs = numpy.append(self.afs, next_af)  # new AF
            self.add_af()  # tell subclass to update internal variables
            # FIXME: S defined in subclass!
            # self.S = numpy.append(self.S, 0) # new S
            self.dev = numpy.append(self.dev, 0)  # new dev
            self.af_idx = next_idx
            self.af_decade = next_decade
        else:
            pass
            # print "no new AF "

    def add_frequency(self, f):
        """ add new frequency point, in units of Hz """
        if not self.x:  # empty sequence
            self.add_phase(0)  # initialize
        self.add_phase(self.x[-1] + f)  # integration

    def taus(self):
        """ return taus, in unit of seconds """
        return self.tau0*numpy.array(self.afs)

    def add_af(self):
        pass  # define in subclass!

    def devs(self):
        """ return deviation """
        return self.dev


class oadev_realtime(dev_realtime):
    """ Overlapping Allan deviation in real-time from a stream of
    phase/frequency samples.

    Reference [Dobrogowski2007]_
    """
    def __init__(self, afs=[1], tau0=1.0, auto_afs=False, pts_per_decade=4):
        super(oadev_realtime, self).__init__(afs=afs, tau0=tau0,
                                             auto_afs=auto_afs,
                                             pts_per_decade=pts_per_decade)
        self.S = numpy.zeros(len(afs))      # sum-of-squares

    def add_phase(self, xnew):
        """ add new phase point, in units of seconds """
        self.x.append(xnew)
        for idx, af in enumerate(self.afs):
            if len(self.x) >= (2*af+1):
                self.update_S(idx)
        if self.auto_afs:
            self.update_af()

    def update_S(self, idx):
        """ update S, sum-of-squares """
        af = self.afs[idx]
        i = len(self.x)-1  # last pt
        S_new = pow(self.x[i] - 2*self.x[i-af] + self.x[i-2*af], 2)
        self.S[idx] = self.S[idx] + S_new
        self.dev[idx] = numpy.sqrt(
            (1.0/(2*pow(af*self.tau0, 2)*(i+1-2*af))) * self.S[idx])

    def add_af(self):
        self.S = numpy.append(self.S, 0)


class ohdev_realtime(dev_realtime):
    """ Overlapping Hadamard deviation in real-time from a stream of
    phase/frequency samples.

    Reference [Dobrogowski2007]_
    """
    def __init__(self, afs=[1], tau0=1.0, auto_afs=False, pts_per_decade=4):
        super(ohdev_realtime, self).__init__(afs=afs, tau0=tau0,
                                             auto_afs=auto_afs,
                                             pts_per_decade=pts_per_decade)
        self.S = numpy.zeros(len(afs))      # sum-of-squares

    def add_af(self):
        self.S = numpy.append(self.S, 0)

    def add_phase(self, xnew):
        """ add new phase point """
        self.x.append(xnew)
        for idx, af in enumerate(self.afs):
            if len(self.x) > 3*af:
                self.update_S(idx)
        if self.auto_afs:
            self.update_af()

    def update_S(self, idx):
        """ update S, sum-of-squares """
        af = self.afs[idx]
        i = len(self.x)-1  # last pt
        # print i,self.x
        S_new = pow(self.x[i] -
                    3*self.x[i-af] +
                    3*self.x[i-2*af] -
                    self.x[i-3*af], 2)
        self.S[idx] = self.S[idx] + S_new
        self.dev[idx] = numpy.sqrt(
            (1.0/(6.0*pow(af*self.tau0, 2)*(i+1.0-3*af))) * self.S[idx])


class tdev_realtime(dev_realtime):
    """ Time deviation and Modified Allan deviation in real-time from a stream
    of phase/frequency samples.

    Reference [Dobrogowski2007]_
    """
    def __init__(self, afs=[1], tau0=1.0, auto_afs=False, pts_per_decade=4):
        super(tdev_realtime, self).__init__(afs=afs, tau0=tau0,
                                            auto_afs=auto_afs,
                                            pts_per_decade=pts_per_decade)
        self.S = numpy.zeros(len(afs))      # sum-of-squares
        self.So = numpy.zeros(len(afs))      # overall sum-of-squares

    def add_phase(self, xnew):
        """ add new phase point """
        self.x.append(xnew)
        for idx, af in enumerate(self.afs):
            if len(self.x) >= 3*af+1:  # 3n+1 samples measured
                self.update_S(idx)
            elif len(self.x) >= 2*af+1:  # 2n+1 samples measured
                self.update_S3n(idx)

        if self.auto_afs:
            self.update_af()

    def add_af(self):
        self.S = numpy.append(self.S, 0)
        self.So = numpy.append(self.So, 0)

    def update_S3n(self, idx):
        """ eqn (13) of paper """
        af = self.afs[idx]
        j = len(self.x)-1  # last pt
        self.S[idx] = self.S[idx] + self.x[j] - 2*self.x[j-af] + self.x[j-2*af]
        if len(self.x) == 3*af:
            # last call to this fctn
            self.So[idx] = pow(self.S[idx], 2)
            self.update_dev(idx)

    def update_dev(self, idx):
        # Eqn (14)
        num_pts = len(self.x)
        af = self.afs[idx]
        self.dev[idx] = numpy.sqrt(
            (1.0/6.0)*(1.0/(num_pts-3*af+1.0))*(1.0/pow(af, 2))*(self.So[idx]))

    def update_S(self, idx):
        """ update S, sum-of-squares """
        af = self.afs[idx]
        assert(len(self.x) >= 3*af+1)
        i = len(self.x)-1  # last pt
        # Eqn (12)
        S_new = (-1*self.x[i-3*af] + 3*self.x[i-2*af] -
                 3*self.x[i-af] + self.x[i])

        self.S[idx] = self.S[idx] + S_new
        # Eqn (11)
        self.So[idx] = self.So[idx] + pow(self.S[idx], 2)
        # ??? S_(i-1) in paper for TDEV-sqrt?
        self.update_dev(idx)

    def mdev(self):
        """ scale tdev to output mdev """
        mdev = self.dev.copy()
        for idx, af in enumerate(self.afs):
            mdev[idx] = mdev[idx]*numpy.sqrt(3)/(af*self.tau0)
        return mdev

# end of file realtime.py

1	"""
2	Real-time statistics
3	====================
4
5	Author: Anders Wallin (anders.e.e.wallin "at" gmail.com)
6
7	Statistics are computed 'on-the-fly' from a stream of
8	phase/frequency samples
9	Initial version 2019 July, Anders E. E. Wallin.
10
11	This file is part of allantools, see https://github.com/aewallin/allantools
12
13	Version history
14	---------------
15
16	2019.07
17	- Initial release
18
19	License
20	-------
21	This program is free software: you can redistribute it and/or modify
22	it under the terms of the GNU Lesser General Public License as published by
23	the Free Software Foundation, either version 3 of the License, or
24	(at your option) any later version.
25	This program is distributed in the hope that it will be useful,
26	but WITHOUT ANY WARRANTY; without even the implied warranty of
27	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
28	GNU Lesser General Public License for more details.
29	You should have received a copy of the GNU Lesser General Public License
30	along with this program. If not, see <http://www.gnu.org/licenses/>.
31	"""
32
33	import numpy	2✔
34
35
36	class dev_realtime(object):	2✔
37	""" Base-class for real-time statistics """
38	def __init__(self, afs=[1], tau0=1.0, auto_afs=False, pts_per_decade=4):	1✔
39	self.x = [] # phase time-series	1✔
40	self.afs = afs # averaging factor, tau = af*tau0	1✔
41	self.auto_afs = auto_afs	1✔
42	# logscpace will fail at >=6 (?), need to remove duplicates?
43	self.af_taus = numpy.logspace(0, 1, pts_per_decade+1)[:-1]	1✔
44	# logspace index, to keep track of afs in auto-af mode
45	self.af_idx = 0	1✔
46	self.af_decade = 0 # logspace decade	1✔
47	if auto_afs:	1✔
48	self.afs = numpy.array([1])	1✔
49	self.dev = numpy.zeros(len(afs)) # resulting xDEV	1✔
50	self.tau0 = tau0 # time-interval between points	1✔
51
52	def update_af(self):	1✔
53	""" used in auto-AF mode,
54	- check if we can add another AF
55	- if yes, add it.
56	"""
57	next_idx = self.af_idx+1	3✔
58	next_decade = self.af_decade	3✔
59	if next_idx == len(self.af_taus):	3✔
60	next_idx = 0	3✔
61	next_decade = self.af_decade + 1	3✔
62
63	# next possible AF:
64	next_af = int(numpy.round(	3✔
65	pow(10.0, next_decade) * self.af_taus[next_idx]))	3✔
66	if len(self.x) >= (2*next_af+1): # can compute next AF	3✔
67	self.afs = numpy.append(self.afs, next_af) # new AF	3✔
68	self.add_af() # tell subclass to update internal variables	3✔
69	# FIXME: S defined in subclass!
70	# self.S = numpy.append(self.S, 0) # new S
71	self.dev = numpy.append(self.dev, 0) # new dev	3✔
72	self.af_idx = next_idx	3✔
73	self.af_decade = next_decade	3✔
74	else:	×
75	pass	3✔
76	# print "no new AF "
77
78	def add_frequency(self, f):	3✔
79	""" add new frequency point, in units of Hz """
80	if not self.x: # empty sequence
81	self.add_phase(0) # initialize
82	self.add_phase(self.x[-1] + f) # integration
83
84	def taus(self):
85	""" return taus, in unit of seconds """
86	return self.tau0*numpy.array(self.afs)	3✔
87
88	def add_af(self):	3✔
89	pass # define in subclass!	×
90
91	def devs(self):	3✔
92	""" return deviation """
93	return self.dev
94
95
96	class oadev_realtime(dev_realtime):
97	""" Overlapping Allan deviation in real-time from a stream of
98	phase/frequency samples.
99
100	Reference [Dobrogowski2007]_
101	"""
102	def __init__(self, afs=[1], tau0=1.0, auto_afs=False, pts_per_decade=4):	3✔
103	super(oadev_realtime, self).__init__(afs=afs, tau0=tau0,	3✔
104	auto_afs=auto_afs,	3✔
105	pts_per_decade=pts_per_decade)	3✔
106	self.S = numpy.zeros(len(afs)) # sum-of-squares	3✔
107
108	def add_phase(self, xnew):	3✔
109	""" add new phase point, in units of seconds """
110	self.x.append(xnew)
111	for idx, af in enumerate(self.afs):
112	if len(self.x) >= (2*af+1):
113	self.update_S(idx)
114	if self.auto_afs:
115	self.update_af()
116
117	def update_S(self, idx):
118	""" update S, sum-of-squares """
119	af = self.afs[idx]	3✔
120	i = len(self.x)-1 # last pt	3✔
121	S_new = pow(self.x[i] - 2self.x[i-af] + self.x[i-2af], 2)	3✔
122	self.S[idx] = self.S[idx] + S_new	3✔
123	self.dev[idx] = numpy.sqrt(	3✔
124	(1.0/(2pow(afself.tau0, 2)(i+1-2af))) * self.S[idx])	3✔
125
126	def add_af(self):	3✔
127	self.S = numpy.append(self.S, 0)	3✔
128
129
130	class ohdev_realtime(dev_realtime):	3✔
131	""" Overlapping Hadamard deviation in real-time from a stream of
132	phase/frequency samples.
133
134	Reference [Dobrogowski2007]_
135	"""
136	def __init__(self, afs=[1], tau0=1.0, auto_afs=False, pts_per_decade=4):	3✔
137	super(ohdev_realtime, self).__init__(afs=afs, tau0=tau0,	3✔
138	auto_afs=auto_afs,	3✔
139	pts_per_decade=pts_per_decade)	3✔
140	self.S = numpy.zeros(len(afs)) # sum-of-squares	3✔
141
142	def add_af(self):	3✔
143	self.S = numpy.append(self.S, 0)	3✔
144
145	def add_phase(self, xnew):	3✔
146	""" add new phase point """
147	self.x.append(xnew)
148	for idx, af in enumerate(self.afs):
149	if len(self.x) > 3*af:
150	self.update_S(idx)
151	if self.auto_afs:
152	self.update_af()
153
154	def update_S(self, idx):
155	""" update S, sum-of-squares """
156	af = self.afs[idx]	3✔
157	i = len(self.x)-1 # last pt	3✔
158	# print i,self.x
159	S_new = pow(self.x[i] -	3✔
160	3*self.x[i-af] +	2✔
161	3self.x[i-2af] -	3✔
162	self.x[i-3*af], 2)	3✔
163	self.S[idx] = self.S[idx] + S_new	3✔
164	self.dev[idx] = numpy.sqrt(	3✔
165	(1.0/(6.0pow(afself.tau0, 2)(i+1.0-3af))) * self.S[idx])	3✔
166
167
168	class tdev_realtime(dev_realtime):	3✔
169	""" Time deviation and Modified Allan deviation in real-time from a stream
170	of phase/frequency samples.
171
172	Reference [Dobrogowski2007]_
173	"""
174	def __init__(self, afs=[1], tau0=1.0, auto_afs=False, pts_per_decade=4):	3✔
175	super(tdev_realtime, self).__init__(afs=afs, tau0=tau0,	3✔
176	auto_afs=auto_afs,	3✔
177	pts_per_decade=pts_per_decade)	3✔
178	self.S = numpy.zeros(len(afs)) # sum-of-squares	3✔
179	self.So = numpy.zeros(len(afs)) # overall sum-of-squares	3✔
180
181	def add_phase(self, xnew):	3✔
182	""" add new phase point """
183	self.x.append(xnew)
184	for idx, af in enumerate(self.afs):
185	if len(self.x) >= 3*af+1: # 3n+1 samples measured
186	self.update_S(idx)
187	elif len(self.x) >= 2*af+1: # 2n+1 samples measured
188	self.update_S3n(idx)
189
190	if self.auto_afs:
191	self.update_af()
192
193	def add_af(self):
194	self.S = numpy.append(self.S, 0)
195	self.So = numpy.append(self.So, 0)
196
197	def update_S3n(self, idx):
198	""" eqn (13) of paper """
199	af = self.afs[idx]	3✔
200	j = len(self.x)-1 # last pt	3✔
201	self.S[idx] = self.S[idx] + self.x[j] - 2self.x[j-af] + self.x[j-2af]	3✔
202	if len(self.x) == 3*af:	3✔
203	# last call to this fctn
204	self.So[idx] = pow(self.S[idx], 2)	3✔
205	self.update_dev(idx)	3✔
206
207	def update_dev(self, idx):	3✔
208	# Eqn (14)
209	num_pts = len(self.x)	3✔
210	af = self.afs[idx]	3✔
211	self.dev[idx] = numpy.sqrt(	3✔
212	(1.0/6.0)(1.0/(num_pts-3af+1.0))(1.0/pow(af, 2))(self.So[idx]))	3✔
213
214	def update_S(self, idx):	3✔
215	""" update S, sum-of-squares """
216	af = self.afs[idx]
217	assert(len(self.x) >= 3*af+1)
218	i = len(self.x)-1 # last pt
219	# Eqn (12)
220	S_new = (-1self.x[i-3af] + 3self.x[i-2af] -
221	3*self.x[i-af] + self.x[i])
222
223	self.S[idx] = self.S[idx] + S_new
224	# Eqn (11)
225	self.So[idx] = self.So[idx] + pow(self.S[idx], 2)
226	# ??? S_(i-1) in paper for TDEV-sqrt?
227	self.update_dev(idx)
228
229	def mdev(self):
230	""" scale tdev to output mdev """
231	mdev = self.dev.copy()	3✔
232	for idx, af in enumerate(self.afs):	3✔
233	mdev[idx] = mdev[idx]numpy.sqrt(3)/(afself.tau0)	3✔
234	return mdev	3✔
235
236	# end of file realtime.py

aewallin / allantools / 6997134815

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous