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

83.33

/allantools/noise.py

"""
Allan deviation tools, Functions for generating noise
=====================================================

See: http://en.wikipedia.org/wiki/Colors_of_noise

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

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

**2019.07**
- Homogeneized output types and parameter names
- PEP8 + docstrings update

**v1.00**
- initial version Anders Wallin, 2014 January


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
   GGNU 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 math
import numpy
import scipy.signal  # for welch PSD


def numpy_psd(x, f_sample=1.0):
    """ calculate power spectral density of input signal x
        x = signal
        f_sample = sampling frequency in Hz. i.e. 1/fs is the time-interval
             in seconds between datapoints
        scale fft so that output corresponds to 1-sided PSD
        output has units of [X^2/Hz] where X is the unit of x
    """
    psd_of_x = ((2.0 / (float(len(x)) * f_sample))
                * numpy.abs(numpy.fft.rfft(x))**2)
    f_axis = numpy.linspace(0, f_sample/2.0, len(psd_of_x))  # frequency axis
    return f_axis, psd_of_x


def scipy_psd(x, f_sample=1.0, nr_segments=4):
    """ PSD routine from scipy
        we can compare our own numpy result against this one
    """
    f_axis, psd_of_x = scipy.signal.welch(x,
                                          f_sample,
                                          nperseg=len(x)/nr_segments)
    return f_axis, psd_of_x


def white(num_points=1024, b0=1.0, fs=1.0):
    """ White noise generator

        Generate time series with white noise that has constant PSD = b0,
        up to the nyquist frequency fs/2.

        The PSD is at 'height' b0 and extends from 0 Hz up to the nyquist
        frequency fs/2 (prefactor math.sqrt(b0*fs/2.0))

        Parameters
        ----------
        num_points: int, optional
            number of samples
        b0: float, optional
            desired power-spectral density in [X^2/Hz] where X is the unit of x
        fs: float, optional
            sampling frequency, i.e. 1/fs is the time-interval between
            datapoints

        Returns
        -------
        White noise sample: numpy.array
    """
    return math.sqrt(b0*fs/2.0)*numpy.random.randn(num_points)


def brown(num_points=1024, b_minus2=1.0, fs=1.0):
    """ Brownian or random walk (diffusion) noise with 1/f^2 PSD

        Not really a color... rather Brownian or random-walk.
        Obtained by integrating white-noise.

        Parameters
        ----------
        num_points: int, optional
            number of samples
        b_minus2: float, optional
            desired power-spectral density is b2*f^-2
        fs: float, optional
            sampling frequency, i.e. 1/fs is the time-interval between
            datapoints

        Returns
        -------
        Random walk sample: numpy.array
    """
    return (1.0/float(fs))*numpy.cumsum(
        white(num_points,
              b0=b_minus2*(4.0*math.pi*math.pi),
              fs=fs))


def violet(num_points=1024, b2=1, fs=1):
    """ Violet noise with f^2 PSD

        Obtained by differentiating white noise

        Parameters
        ----------
        num_points: int, optional
            number of samples
        b2: float, optional
            desired power-spectral density is b2*f^2
        fs: float, optional
            sampling frequency, i.e. 1/fs is the time-interval between
            datapoints

        Returns
        -------
        Violet noise sample: numpy.array
    """
    # diff() reduces number of points by one.
    return (float(fs))*numpy.diff(
        white(num_points+1, b0=b2/(2.0*math.pi)**2, fs=fs))


def pink(num_points=1024, depth=80):
    """ Pink noise (approximation) with 1/f PSD

        Fills a sample with results from a pink noise generator
        from http://pydoc.net/Python/lmj.sound/0.1.1/lmj.sound.noise/,
        based on the Voss-McCartney algorithm, discussion and code examples at
        http://www.firstpr.com.au/dsp/pink-noise/

        Parameters
        ----------
        num_points: int, optional
            number of samples
        depth: int, optional
            number of iteration for each point. High numbers are slower but
            generates a more correct spectrum on low-frequencies end.

        Returns
        -------
        Pink noise sample: numpy.array
    """
    # FIXME: couldn't we implement here the normalization as for the other
    # noise types using a noise power law coefficient b_minus1?
    a = []
    s = iterpink(depth)
    for n in range(num_points):  # FIXME: num_points is unused here.
        a.append(next(s))
    return numpy.array(a)


def iterpink(depth=20):
    """Generate a sequence of samples of pink noise.

    pink noise generator
    from http://pydoc.net/Python/lmj.sound/0.1.1/lmj.sound.noise/

    Based on the Voss-McCartney algorithm, discussion and code examples at
    http://www.firstpr.com.au/dsp/pink-noise/

    depth: Use this many samples of white noise to calculate the output. A
      higher  number is slower to run, but renders low frequencies with more
      correct power spectra.

    Generates a never-ending sequence of floating-point values. Any continuous
    set of these samples will tend to have a 1/f power spectrum.
    """
    values = numpy.random.randn(depth)
    smooth = numpy.random.randn(depth)
    source = numpy.random.randn(depth)
    sumvals = values.sum()
    i = 0
    while True:
        yield sumvals + smooth[i]

        # advance the index by 1. if the index wraps, generate noise to use in
        # the calculations, but do not update any of the pink noise values.
        i += 1
        if i == depth:
            i = 0
            smooth = numpy.random.randn(depth)
            source = numpy.random.randn(depth)
            continue

        # count trailing zeros in i
        c = 0
        while not (i >> c) & 1:
            c += 1

        # replace value c with a new source element
        sumvals += source[i] - values[c]
        values[c] = source[i]

1	"""
2	Allan deviation tools, Functions for generating noise
3	=====================================================
4
5	See: http://en.wikipedia.org/wiki/Colors_of_noise
6
7	Author: Anders Wallin (anders.e.e.wallin "at" gmail.com)
8
9	Version history
10	---------------
11
12	2019.07
13	- Homogeneized output types and parameter names
14	- PEP8 + docstrings update
15
16	v1.00
17	- initial version Anders Wallin, 2014 January
18
19
20	This program is free software: you can redistribute it and/or modify
21	it under the terms of the GNU Lesser General Public License as published by
22	the Free Software Foundation, either version 3 of the License, or
23	(at your option) any later version.
24
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	GGNU Lesser General Public License for more details.
29
30	You should have received a copy of the GNU Lesser General Public License
31	along with this program. If not, see <http://www.gnu.org/licenses/>.
32	"""
33
34	import math	2✔
35	import numpy	2✔
36	import scipy.signal # for welch PSD	2✔
37
38
39	def numpy_psd(x, f_sample=1.0):	2✔
40	""" calculate power spectral density of input signal x
41	x = signal
42	f_sample = sampling frequency in Hz. i.e. 1/fs is the time-interval
43	in seconds between datapoints
44	scale fft so that output corresponds to 1-sided PSD
45	output has units of [X^2/Hz] where X is the unit of x
46	"""
47	psd_of_x = ((2.0 / (float(len(x)) * f_sample))	×
48	* numpy.abs(numpy.fft.rfft(x))**2)	×
49	f_axis = numpy.linspace(0, f_sample/2.0, len(psd_of_x)) # frequency axis	×
50	return f_axis, psd_of_x	×
51
52
53	def scipy_psd(x, f_sample=1.0, nr_segments=4):	3✔
54	""" PSD routine from scipy
55	we can compare our own numpy result against this one
56	"""
57	f_axis, psd_of_x = scipy.signal.welch(x,	×
58	f_sample,	×
59	nperseg=len(x)/nr_segments)	×
60	return f_axis, psd_of_x	×
61
62
63	def white(num_points=1024, b0=1.0, fs=1.0):	3✔
64	""" White noise generator
65
66	Generate time series with white noise that has constant PSD = b0,
67	up to the nyquist frequency fs/2.
68
69	The PSD is at 'height' b0 and extends from 0 Hz up to the nyquist
70	frequency fs/2 (prefactor math.sqrt(b0*fs/2.0))
71
72	Parameters
73	----------
74	num_points: int, optional
75	number of samples
76	b0: float, optional
77	desired power-spectral density in [X^2/Hz] where X is the unit of x
78	fs: float, optional
79	sampling frequency, i.e. 1/fs is the time-interval between
80	datapoints
81
82	Returns
83	-------
84	White noise sample: numpy.array
85	"""
86	return math.sqrt(b0fs/2.0)numpy.random.randn(num_points)	3✔
87
88
89	def brown(num_points=1024, b_minus2=1.0, fs=1.0):	3✔
90	""" Brownian or random walk (diffusion) noise with 1/f^2 PSD
91
92	Not really a color... rather Brownian or random-walk.
93	Obtained by integrating white-noise.
94
95	Parameters
96	----------
97	num_points: int, optional
98	number of samples
99	b_minus2: float, optional
100	desired power-spectral density is b2*f^-2
101	fs: float, optional
102	sampling frequency, i.e. 1/fs is the time-interval between
103	datapoints
104
105	Returns
106	-------
107	Random walk sample: numpy.array
108	"""
109	return (1.0/float(fs))*numpy.cumsum(	3✔
110	white(num_points,	3✔
111	b0=b_minus2(4.0math.pi*math.pi),	3✔
112	fs=fs))	3✔
113
114
115	def violet(num_points=1024, b2=1, fs=1):	3✔
116	""" Violet noise with f^2 PSD
117
118	Obtained by differentiating white noise
119
120	Parameters
121	----------
122	num_points: int, optional
123	number of samples
124	b2: float, optional
125	desired power-spectral density is b2*f^2
126	fs: float, optional
127	sampling frequency, i.e. 1/fs is the time-interval between
128	datapoints
129
130	Returns
131	-------
132	Violet noise sample: numpy.array
133	"""
134	# diff() reduces number of points by one.
135	return (float(fs))*numpy.diff(	3✔
136	white(num_points+1, b0=b2/(2.0math.pi)*2, fs=fs))	3✔
137
138
139	def pink(num_points=1024, depth=80):	3✔
140	""" Pink noise (approximation) with 1/f PSD
141
142	Fills a sample with results from a pink noise generator
143	from http://pydoc.net/Python/lmj.sound/0.1.1/lmj.sound.noise/,
144	based on the Voss-McCartney algorithm, discussion and code examples at
145	http://www.firstpr.com.au/dsp/pink-noise/
146
147	Parameters
148	----------
149	num_points: int, optional
150	number of samples
151	depth: int, optional
152	number of iteration for each point. High numbers are slower but
153	generates a more correct spectrum on low-frequencies end.
154
155	Returns
156	-------
157	Pink noise sample: numpy.array
158	"""
159	# FIXME: couldn't we implement here the normalization as for the other
160	# noise types using a noise power law coefficient b_minus1?
161	a = []	3✔
162	s = iterpink(depth)	3✔
163	for n in range(num_points): # FIXME: num_points is unused here.	3✔
164	a.append(next(s))	3✔
165	return numpy.array(a)	3✔
166
167
168	def iterpink(depth=20):	3✔
169	"""Generate a sequence of samples of pink noise.
170
171	pink noise generator
172	from http://pydoc.net/Python/lmj.sound/0.1.1/lmj.sound.noise/
173
174	Based on the Voss-McCartney algorithm, discussion and code examples at
175	http://www.firstpr.com.au/dsp/pink-noise/
176
177	depth: Use this many samples of white noise to calculate the output. A
178	higher number is slower to run, but renders low frequencies with more
179	correct power spectra.
180
181	Generates a never-ending sequence of floating-point values. Any continuous
182	set of these samples will tend to have a 1/f power spectrum.
183	"""
184	values = numpy.random.randn(depth)	3✔
185	smooth = numpy.random.randn(depth)	3✔
186	source = numpy.random.randn(depth)	3✔
187	sumvals = values.sum()	3✔
188	i = 0	3✔
189	while True:	1✔
190	yield sumvals + smooth[i]	3✔
191
192	# advance the index by 1. if the index wraps, generate noise to use in
193	# the calculations, but do not update any of the pink noise values.
194	i += 1	3✔
195	if i == depth:	3✔
196	i = 0	3✔
197	smooth = numpy.random.randn(depth)	3✔
198	source = numpy.random.randn(depth)	3✔
199	continue	3✔
200
201	# count trailing zeros in i
202	c = 0	3✔
203	while not (i >> c) & 1:	3✔
204	c += 1	3✔
205
206	# replace value c with a new source element
207	sumvals += source[i] - values[c]	3✔
208	values[c] = source[i]	3✔

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