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
|
import scipy.special
from scipy_base import *
from scipy_base.fastumath import sqrt, exp, greater, equal, cos, add, sin
from spline import * # C-modules
gamma = scipy.special.gamma
def factorial(n):
return gamma(n+1)
def spline_filter(Iin, lmbda=5.0):
"""Smoothing spline (cubic) filtering of a rank-2 array.
Filter an input data set, Iin, using a (cubic) smoothing spline of
fall-off lmbda.
"""
intype = Iin.typecode()
hcol = sarray([1.0,4.0,1.0],'f')/6.0
if intype in ['F','D']:
Iin = Iin.astype('F')
ckr = cspline2d(Iin.real,lmbda)
cki = cspline2d(Iin.imag,lmbda)
outr = sepfir2d(ckr,hcol,hcol)
outi = sepfir2d(cki,hcol,hcol)
out = (outr + 1j*outi).astype(intype)
elif intype in ['f','d']:
ckr = cspline2d(Iin,lmbda)
out = sepfir2d(ckr, hcol, hcol)
out = out.astype(intype)
else:
raise TypeError;
return out
def _bspline(x,n):
"""bspline(x,n) -> y: B-spline basis function of order n.
"""
jlist = arange(n+2)
val = 0.0
baseval = x + (n+1)/2.0
for j in jlist:
xval = baseval - j
if xval >=0:
if j % 2 == 0:
fact = 1
else:
fact = -1
term = fact * (n+1) * xval**n / factorial(j) / factorial(n+1-j)
val = val + term
return val
bspline = vectorize(_bspline)
def gauss_spline(x,n):
"""Gaussian approximation to B-spline basis function of order n.
"""
signsq = (n+1) / 12.0
return 1/sqrt(2*pi*signsq) * exp(-x**2 / 2 / signsq)
def cubic(x):
ax = abs(x)
f1 = 2.0/3 - 1.0/2*ax**2 * (2-ax)
f2 = 1.0/6*(2-ax)**3
f3 = where(greater(ax,1),f2,f1)
return where(greater(ax,2),0,f3)
def quintic(n):
n = asarray(n)
an = abs(n)
f1 = where(equal(an,0),66.0/120.0,0)
f2 = where(equal(an,1),26.0/120.0,0)
f3 = where(equal(an,2),1.0/120.0,0)
return f1 + f2 + f3
def quadratic(x):
ax = abs(x)
f1 = 0.75-ax**2
f2 = (ax-1.5)**2 / 2.0
f3 = where(greater(abs(x),0.5),f2,f1)
return where(greater(abs(x),1.5),0,f3)
def c0_P(order):
# values taken from Unser, et.al. 1993 IEEE
if order == 0:
c0 = 1
P = array([1])
elif order == 1:
c0 = 1
P = array([0,1])
elif order == 2:
c0 = 8
P = array([1,6,1])
elif order == 3:
c0 = 6
P = array([1,4,1])
elif order == 4:
c0 = 384
P = array([1,76,230,76,1])
elif order == 5:
c0 = 120
P = array([1,26,66,26,1])
elif order == 6:
c0 = 46080
P = array([1,722,10543,23548, 10543, 722, 1])
elif order == 7:
c0 = 5040
P = array([1,120,1191,2416,1191, 120, 1])
else:
raise ValueError, "Unknown order."
def _coeff_smooth(lam):
xi = 1 - 96*lam + 24*lam * sqrt(3 + 144*lam)
omeg = arctan2(sqrt(144*lam-1),sqrt(xi))
rho = (24*lam - 1 - sqrt(xi)) / (24*lam)
rho = rho * sqrt((48*lam + 24*lam * sqrt(3+144*lam))/xi)
return rho,omeg
def _cubic_smooth_coeff(signal,lamb):
rho, omega = _coeff_smooth(lamb)
cs = 1-2*rho*cos(omega) + rho*rho
K = len(signal)
yp = zeros((K,),signal.typecode())
k = arange(K)
yp[0] = hc(0,cs,rho,omega)*signal[0] + \
add.reduce(hc(k+1,cs,rho,omega)*signal)
yp[1] = hc(0,cs,rho,omega)*signal[0] + \
hc(1,cs,rho,omega)*signal[1] + \
add.reduce(hc(k+2,cs,rho,omega)*signal)
for n in range(2,K):
yp[n] = cs * signal[n] + 2*rho*cos(omega)*yp[n-1] - rho*rho*yp[n-2]
y = zeros((K,),signal.typecode())
y[K-1] = add.reduce((hs(k,cs,rho,omega) + hs(k+1,cs,rho,omega))*signal[::-1])
y[K-2] = add.reduce((hs(k-1,cs,rho,omega) + hs(k+2,cs,rho,omega))*signal[::-1])
for n in range(K-3,-1,-1):
y[n] = cs*yp[n] + 2*rho*cos(omega)*y[n+1] - rho*rho*y[n+2]
return y
def _cubic_coeff(signal):
zi = -2 + sqrt(3)
K = len(signal)
yplus = zeros((K,),signal.typecode())
powers = zi**arange(K)
yplus[0] = signal[0] + zi*add.reduce(powers*signal)
for k in range(1,K):
yplus[k] = signal[k] + zi*yplus[k-1]
output = zeros((K,),signal.typecode())
output[K-1] = zi / (zi-1)*yplus[K-1]
for k in range(K-2,-1,-1):
output[k] = zi*(output[k+1]-yplus[k])
return output*6.0
def _quadratic_coeff(signal):
zi = -3 + 2*sqrt(2.0)
K = len(signal)
yplus = zeros((K,),signal.typecode())
powers = zi**arange(K)
yplus[0] = signal[0] + zi*add.reduce(powers*signal)
for k in range(1,K):
yplus[k] = signal[k] + zi*yplus[k-1]
output = zeros((K,),signal.typecode())
output[K-1] = zi / (zi-1)*yplus[K-1]
for k in range(K-2,-1,-1):
output[k] = zi*(output[k+1]-yplus[k])
return output*8.0
def cspline1d(signal,lamb=0.0):
"""Compute cubic spline coefficients for rank-1 array.
Description:
Find the cubic spline coefficients for a 1-D signal assuming
mirror-symmetric boundary conditions. To obtain the signal back from
the spline representation mirror-symmetric-convolve these coefficients
with a length 3 FIR window [1.0, 4.0, 1.0]/ 6.0 .
Inputs:
signal -- a rank-1 array representing samples of a signal.
lamb -- smoothing coefficient (default = 0.0)
Output:
c -- cubic spline coefficients.
"""
if lamb != 0.0:
return _cubic_smooth_coeff(signal,lamb)
else:
return _cubic_coeff(signal)
def qspline1d(signal,lamb=0.0):
"""Compute quadratic spline coefficients for rank-1 array.
Description:
Find the quadratic spline coefficients for a 1-D signal assuming
mirror-symmetric boundary conditions. To obtain the signal back from
the spline representation mirror-symmetric-convolve these coefficients
with a length 3 FIR window [1.0, 6.0, 1.0]/ 8.0 .
Inputs:
signal -- a rank-1 array representing samples of a signal.
lamb -- smoothing coefficient (must be zero for now.)
Output:
c -- cubic spline coefficients.
"""
if lamb != 0.0:
raise ValueError, "Smoothing quadratic splines not supported yet."
else:
return _quadratic_coeff(signal)
def hc(k,cs,rho,omega):
return cs / sin(omega) * (rho**k)*sin(omega*(k+1))*(greater(k,-1))
def hs(k,cs,rho,omega):
c0 = cs*cs * (1 + rho*rho) / (1 - rho*rho) / (1-2*rho*rho*cos(2*omega) + rho**4)
gamma = (1-rho*rho) / (1+rho*rho) / tan(omega)
ak = abs(k)
return c0 * rho**ak * (cos(omega*ak) + gamma*sin(omega*ak))
|
