import sys
import csv
import matplotlib.pyplot as plt
from math import sqrt, exp, pi
import matplotlib.pyplot as plt
"""
Smoothing by convolution with Gauss function
"""
#======================
# parameters
#======================
csvfile = 'dos.csv'
outfile = 'convolution.csv'
width = 0.2
if __name__ == "__main__":
argv = sys.argv
n = len(argv)
if n >= 2:
width = float(argv[1])
[ドキュメント]
def convolution(x, y, width):
ndata = len(x)
dx = x[1] - x[0]
xrange = x[ndata-1] - x[0]
# integration range, converted to number of the list index
di = int( (width * 5.0) / dx + 1.1 )
# the coefficient of Gauss function
coeff = 1.0 / sqrt(pi)/ width
# deconvoluted data
ys = [0.0]*ndata
for j in range(0, ndata):
x0 = x[j];
y0 = y[j];
for k in range(-di, di+1):
if j+k < 0 or j+k >= ndata:
continue
dvx = dx * k / width
f = dx * coeff * exp(-dvx*dvx)
ys[j+k] += y0 * f
return ys;
[ドキュメント]
def savecsv(outfile, x, y, ys):
try:
print("Write to [{}]".format(outfile))
f = open(outfile, 'w')
except:
# except IOError:
print("Error: Can not write to [{}]".format(outfile))
else:
fout = csv.writer(f, lineterminator='\n')
fout.writerow(('x', 'y(raw)', 'y(smooth)'))
# fout.writerows(data)
for i in range(0, len(x)):
fout.writerow((x[i], y[i], ys[i]))
f.close()
[ドキュメント]
def main():
global csvfile
global outfile
global width
print("Read data from [{}]".format(csvfile))
x = []
y = []
with open(csvfile) as f:
fin = csv.reader(f)
next(fin)
c = 0
for row in fin:
if row[0] == '':
break
if c == 0:
label = row
print("{}\t{}".format(label[0], label[1]))
else:
print("{}\t{}".format(row[0], row[1]))
x.append(float(row[0]))
y.append(float(row[1]))
c += 1
ndata = len(x)
print("")
print("Smoothing by convolution with Gauss function of w = {}".format(width))
ys = convolution(x, y, width)
for i in range(0, ndata):
print("{}\t{}\t{}".format(x[i], y[i], ys[i]))
savecsv(outfile, x, y, ys)
print("")
#=============================
# Plot graphs
#=============================
fig = plt.figure()
ax1 = fig.add_subplot(1, 1, 1)
ax1.plot(x, y, label = 'raw data', linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 0.5)
ax1.plot(x, ys, label = 'convoluted')
ax1.set_xlabel("x")
ax1.set_ylabel("y")
ax1.legend()
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
exit()
if __name__ == '__main__':
main()