-
Notifications
You must be signed in to change notification settings - Fork 54
/
_synthesizeNTF0.py
287 lines (264 loc) · 9.55 KB
/
_synthesizeNTF0.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
# -*- coding: utf-8 -*-
# _synthesizeNTF0.py
# Module providing the synthesizeNTF function
# Copyright 2013 Giuseppe Venturini
# This file is distributed with python-deltasigma.
#
# python-deltasigma is a 1:1 Python replacement of Richard Schreier's
# MATLAB delta sigma toolbox (aka "delsigma"), upon which it is heavily based.
# The delta sigma toolbox is (c) 2009, Richard Schreier.
#
# python-deltasigma 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
# LICENSE file for the licensing terms.
#
# The following code has been (slightly) modified from pydsm, its original
# copyright notice follows:
#
# Copyright (c) 2012, Sergio Callegari
# All rights reserved.
#
# The py dsm code was ported from the MATLAB Delta Sigma toolbox
# Copyright (c) 2009, Richard Schreier
#
# The three software follow the same license, known as the 3-clause BSD.
# See the LICENSE file for details.
"""
Synthesize a noise transfer function (NTF) for a delta-sigma modulator without
optimizing the result.
"""
import numpy as np
from warnings import warn
from ._evalTF import evalTF
from ._utils import cplxpair, zpk
from ._ds_optzeros import ds_optzeros
from ._config import itn_limit
def synthesizeNTF0(order, osr, opt, H_inf, f0):
"""Synthesize a noise transfer function for a delta-sigma modulator.
::warn This function is not meant to be used directly, instead, set
optimize_NTF to False and call synthesizeNTF(...).
Parameters
----------
order : int,
the order of the modulator
osr : float,
the oversamping ratio
opt : int or list of floats
flag for optimized zeros
* 0 -> not optimized,
* 1 -> optimized,
* 2 -> optimized with at least one zero at band-center,
* [z] -> zero locations in complex form
H_inf : real
max allowed peak value of the NTF.
f0 : real
center frequency for BP modulators, or 0 for LP modulators.
Note that 1 corresponds to the sampling frequency, so that 0.5 is the
maximum value. Value 0 specifies an LP modulator.
Returns
-------
ntf : tuple
noise transfer function in zpk form.
Raises
------
ValueError
'Cannot synthesize NTF zeros' if an empty list is supplied for the zeros.
'Order must be even for a bandpass modulator.' if the order is
incompatible with the modulator type.
'The opt vector must be of length xxx' if opt is used to explicitly
pass the NTF zeros and these are in the wrong number.
Warns
-----
'Unable to achieve specified H_inf ...' if the desired H_inf
cannot be achieved.
'Iteration limit exceeded' if the routine converges too slowly.
Notes
-----
This is actually a wrapper function which calls the appropriate version
of synthesizeNTF, based on the module control flag `optimize_NTF` which
determines whether to use optimization tools.
Parameter H_inf is used to enforce the Lee stability criterion.
See also:
clans() "Closed-loop analysis of noise-shaper." An alternative
method for selecting NTFs based on the 1-norm of the
impulse response of the NTF
synthesizeChebyshevNTF() Select a type-2 highpass Chebyshev NTF.
This function does a better job than synthesizeNTF if osr
or H_inf is low.
"""
# Determine the zeros.
if f0 != 0:
# Bandpass design-- halve the order temporarily.
order = order/2
dw = np.pi/(2*osr)
else:
dw = np.pi/osr
if np.isscalar(opt):
# opt is a number
if opt == 0:
z = np.zeros(order)
else:
z = dw*ds_optzeros(order, opt)
if z.size == 0:
raise ValueError('Cannot synthesize NTF zeros')
if f0 != 0:
# Bandpass design-- shift and replicate the zeros.
order = order*2
z = z + 2*np.pi*f0
z = np.vstack((z, -z)).transpose().flatten()
z = np.exp(1j*z)
else:
z = opt
p = np.zeros(order)
k = 1
fprev = 0
if f0 == 0:
# Lowpass design
HinfLimit = 2**order
# !!! The limit is actually lower for opt=1 and low OSR
if H_inf >= HinfLimit:
warn('Unable to achieve specified H_inf.\n'
'Setting all NTF poles to zero.')
p = np.zeros(order)
else:
x = 0.3**(order-1) # starting guess
for itn in range(1, itn_limit + 1):
me2 = -0.5*(x**(2./order))
w = (2*np.arange(1, order + 1) + 1)*np.pi/order
mb2 = 1 + me2*np.exp(1j*w)
p = mb2 - np.sqrt(mb2**2 - 1)
# Reflect poles to be inside the unit circle
out = abs(p) > 1
p[out] = 1./p[out]
# The following is not exactly what delsig does.
# We do not have an identical cplxpair
p = cplxpair(p)
f = np.real(evalTF(zpk(z, p, k), -1))-H_inf
if itn == 1:
delta_x = -f/100.
else:
delta_x = -f*delta_x/(f - fprev)
xplus = x + delta_x
if xplus > 0:
x = xplus
else:
x = x*0.1
fprev = f
if abs(f) < 1e-10 or abs(delta_x) < 1e-10:
break
if x > 1e6:
warn('Unable to achieve specified Hinf.\n'
'Setting all NTF poles to zero.')
p = np.zeros(order)
break
if itn == itn_limit:
warn('Iteration limit exceeded.')
else:
# Bandpass design
x = 0.3**(order/2-1) # starting guess (not very good for f0~0)
if f0 > 0.25:
z_inf = 1.
else:
z_inf = -1.
c2pif0 = np.cos(2*np.pi*f0)
for itn in range(1, itn_limit+1):
e2 = 0.5*x**(2./order)
w = (2*np.arange(order)+1)*np.pi/order
mb2 = c2pif0 + e2*np.exp(1j*w)
p = mb2 - np.sqrt(mb2**2-1)
# Reflect poles to be inside the unit circle
out = abs(p)>1
p[out] = 1/p[out]
# The following is not exactly what delsig does.
p = cplxpair(p)
f = np.real(evalTF((z, p, k), z_inf)) - H_inf
if itn == 1:
delta_x = -f/100
else:
delta_x = -f*delta_x/(f - fprev)
xplus = x + delta_x
if xplus > 0:
x = xplus
else:
x = x*0.1
fprev = f
if abs(f) < 1e-10 or abs(delta_x) < 1e-10:
break
if x > 1e6:
warn('Unable to achieve specified Hinf.\n'
'Setting all NTF poles to zero.')
p = np.zeros(order)
break
if itn == itn_limit:
warn('Iteration limit exceeded.')
z = cplxpair(z)
return z, p, k
def test_synthesizeNTF0():
"""Test unit for synthesizeNTF0"""
from ._utils import cplxpair
from ._config import optimize_NTF
z, p, k = synthesizeNTF0(order=3, osr=64, opt=0, H_inf=1.5, f0=0.0)
zref = [1., 1., 1.]
pref = [.6694, .7654 + .2793j, .7654 - .2793j]
kref = 1.
assert np.allclose(cplxpair(z), cplxpair(zref), atol=1e-4, rtol=1e-4)
assert np.allclose(cplxpair(p), cplxpair(pref), atol=1e-4, rtol=1e-4)
assert np.allclose(k, kref, atol=1e-4, rtol=1e-4)
# Up next: even order bandpass test
z, p, k = synthesizeNTF0(order=4, osr=32, opt=0, H_inf=1.3, f0=.33)
zref = [-0.4818 + 0.8763j, -0.4818 - 0.8763j, -0.4818 + 0.8763j,
-0.4818 - 0.8763j]
pref = [-0.5125 - 0.7018j, -0.5125 + 0.7018j, -0.3233 - 0.8240j,
-0.3233 + 0.8240j]
kref = 1.
assert np.allclose(cplxpair(z), cplxpair(zref), atol=1e-4, rtol=1e-4)
assert np.allclose(cplxpair(p), cplxpair(pref), atol=1e-4, rtol=1e-4)
assert np.allclose(k, kref, atol=1e-4, rtol=1e-4)
# repeat with zeros optimization
z, p, k = synthesizeNTF0(order=3, osr=64, opt=1, H_inf=1.5, f0=0.0)
zref = [1.0000 + 0.0000j, 0.9993 + 0.0380j, 0.9993 - 0.0380j]
pref = [0.7652 - 0.2795j, 0.7652 + 0.2795j, 0.6692 + 0.0000j]
kref = 1.
assert np.allclose(cplxpair(z), cplxpair(zref), atol=1e-4, rtol=1e-4)
assert np.allclose(cplxpair(p), cplxpair(pref), atol=1e-4, rtol=1e-4)
assert np.allclose(k, kref, atol=1e-4, rtol=1e-4)
# Up next: even order bandpass test
z, p, k = synthesizeNTF0(order=4, osr=32, opt=1, H_inf=1.3, f0=.33)
zref = [-0.4567 + 0.8896j, -0.4567 - 0.8896j, -0.5064 + 0.8623j,
-0.5064 - 0.8623j]
pref = [-0.5125 - 0.7014j, -0.5125 + 0.7014j, -0.3230 - 0.8239j,
-0.3230 + 0.8239j]
kref = 1.
assert np.allclose(cplxpair(z), cplxpair(zref), atol=1e-4, rtol=1e-4)
assert np.allclose(cplxpair(p), cplxpair(pref), atol=1e-4, rtol=1e-4)
assert np.allclose(k, kref, atol=1e-4, rtol=1e-4)
z, p, k = synthesizeNTF0(order=3, osr=64, opt=2, H_inf=1.5, f0=0.0)
zref = [1.0000 + 0.0000j, 0.9993 + 0.0380j, 0.9993 - 0.0380j]
pref = [0.7652 - 0.2795j, 0.7652 + 0.2795j, 0.6692 + 0.0000j]
kref = 1.
assert np.allclose(cplxpair(z), cplxpair(zref), atol=1e-4, rtol=1e-4)
assert np.allclose(cplxpair(p), cplxpair(pref), atol=1e-4, rtol=1e-4)
assert np.allclose(k, kref, atol=1e-4, rtol=1e-4)
# Up next: even order bandpass test
z, p, k = synthesizeNTF0(order=4, osr=32, opt=2, H_inf=1.3, f0=.33)
zref = [-0.4818 + 0.8763j, -0.4818 - 0.8763j, -0.4818 + 0.8763j,
-0.4818 - 0.8763j]
pref = [-0.5125 - 0.7018j, -0.5125 + 0.7018j, -0.3233 - 0.8240j,
-0.3233 + 0.8240j]
kref = 1.
assert np.allclose(cplxpair(z), cplxpair(zref), atol=1e-4, rtol=1e-4)
assert np.allclose(cplxpair(p), cplxpair(pref), atol=1e-4, rtol=1e-4)
assert np.allclose(k, kref, atol=1e-4, rtol=1e-4)
# zeros passed explicitly
opt = [1.0000 + 0.0000j, 0.9986 + 0.06j, 0.9986 - 0.06j,
0.9960 + 0.0892j, 0.9960 - 0.0892j]
z, p, k = synthesizeNTF0(order=5, osr=32, opt=opt, H_inf=1.3, f0=0.0)
zref = [1.0000 + 0.0000j, 0.9986 + 0.06j, 0.9986 - 0.06j,
0.9960 + 0.0892j, 0.9960 - 0.0892j]
pref = [0.8718 - 0.0840j, 0.8718 + 0.0840j, 0.9390 - 0.1475j,
0.9390 + 0.1475j, 0.8491 + 0.0000j]
kref = 1.
assert np.allclose(cplxpair(z), cplxpair(zref), atol=1e-4, rtol=1e-3)
assert np.allclose(cplxpair(p), cplxpair(pref), atol=1e-4, rtol=1e-3)
assert np.allclose(k, kref, atol=1e-4, rtol=1e-4)