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
|
# -*- coding: utf-8 -*-
# We are trying to use fmin to make the deconvolution.
import scipy.optimize as opt
import numpy as np
def test_fmin_convolve(fminfunc, x, h, y, yn, x0):
'''
(*) represents the convolution.
yn is the result of y + noise.
x0 is the initial value of x.
'''
def convolve_func(h):
'''
Calculate the (yn - x(*)h)**2
and make the result the minimum by fmin
'''
return np.sum((yn - np.convolve(x, h))**2)
# call fmin and initialize with x0
h0 = fminfunc(convolve_func, x0)
print(fminfunc.__name__)
print("---------------------")
# relative error between x(*)h0 and y
print("error of y:", np.sum((np.convolve(x, h0)-y)**2)/np.sum(y**2))
print("error of h:", np.sum((h0-h)**2)/np.sum(h**2))
print
def test_n(m, n, nscale):
'''
x, h, y, yn, x0 by random
m is the length of x
n is the length of h
nscale is the intensity of interference.
'''
x = np.random.rand(m)
h = np.random.rand(n)
y = np.convolve(x, h)
yn = y +np.random.rand(len(y)) * nscale
x0 = np.random.rand(n)
test_fmin_convolve(opt.fmin, x, h, y, yn, x0)
test_fmin_convolve(opt.fmin_powell, x, h, y, yn, x0)
test_fmin_convolve(opt.fmin_cg, x, h, y, yn, x0)
test_fmin_convolve(opt.fmin_bfgs, x, h, y, yn, x0)
if __name__ == "__main__":
test_n(200, 20, 0.1)
|
近期评论