104 lines
2.5 KiB
Python
104 lines
2.5 KiB
Python
from scipy import signal
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
import pylab
|
|
from IPython import embed
|
|
from scipy.optimize import curve_fit
|
|
from jar_functions import sin_response
|
|
|
|
def take_second(elem):
|
|
return elem[1]
|
|
|
|
predict = []
|
|
|
|
gain = []
|
|
mgain = []
|
|
phaseshift = []
|
|
mphaseshift = []
|
|
amfreq = []
|
|
amf = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1]
|
|
|
|
currf = None
|
|
idxlist = []
|
|
|
|
data = sorted(np.load('files.npy'), key = take_second)
|
|
|
|
for i, d in enumerate(data):
|
|
dd = list(d)
|
|
jar = np.load('%s.npy' %dd)
|
|
jm = jar - np.mean(jar)
|
|
print(dd)
|
|
|
|
time = np.load('time: %s.npy' %dd)
|
|
|
|
b, a = signal.butter(4, (float(d[1]) / 2) / 10000, 'high', analog=True)
|
|
y = signal.filtfilt(b, a, jm)
|
|
#plt.plot(time, y)
|
|
#plt.plot(time, jar)
|
|
|
|
sinv, sinc = curve_fit(sin_response, time, y, [float(d[1]), 2, 0.5])
|
|
print('frequency, phaseshift, amplitude:', sinv)
|
|
|
|
phaseshift.append(np.sqrt(sinv[1]**2))
|
|
gain.append(np.sqrt(sinv[2]**2))
|
|
amfreq.append(d[1])
|
|
|
|
Rs = []
|
|
for ix, t in enumerate(time):
|
|
R = (jm[ix] - sin_response(t, float(d[1]), np.sqrt(sinv[1]**2), np.sqrt(sinv[2]**2)))**2
|
|
Rs.append(R)
|
|
|
|
sigma = sum(Rs)
|
|
rms = np.sqrt((1/len(time)) * sigma)
|
|
#plt.plot(time, sin_response(time, *sinv), label='fit: f=%f, p=%.2f, A=%.2f' % tuple(sinv))
|
|
|
|
#mean over same amfreqs for phase and gain
|
|
if currf is None or currf == d[1]:
|
|
currf = d[1]
|
|
idxlist.append(i)
|
|
|
|
else: # currf != f
|
|
meanf = [] # lists to make mean of
|
|
meanp = []
|
|
for x in idxlist:
|
|
meanf.append(gain[x])
|
|
meanp.append(phaseshift[x])
|
|
meanedf = np.mean(meanf)
|
|
meanedp = np.mean(meanp)
|
|
mgain.append(meanedf)
|
|
mphaseshift.append(meanedp)
|
|
currf = d[1] # set back for next loop
|
|
idxlist = [i]
|
|
|
|
meanf = []
|
|
meanp = []
|
|
for y in idxlist:
|
|
meanf.append(gain[y])
|
|
meanp.append(phaseshift[y])
|
|
meanedf = np.mean(meanf)
|
|
meanedp = np.mean(meanp)
|
|
mgain.append(meanedf)
|
|
mphaseshift.append(meanedp)
|
|
|
|
for f in amf:
|
|
G = np.max(mgain) / np.sqrt(1 + (2*((np.pi*f*3.14)**2)))
|
|
predict.append(G)
|
|
|
|
|
|
|
|
fig = plt.figure()
|
|
ax = fig.add_subplot(1, 1, 1)
|
|
ax.plot(amf, mgain, 'o')
|
|
ax.plot(amf, predict)
|
|
ax.set_yscale('log')
|
|
ax.set_xscale('log')
|
|
ax.set_title('2018lepto98')
|
|
ax.set_ylabel('gain [Hz/(mV/cm)]')
|
|
ax.set_xlabel('AM-frequency [Hz]')
|
|
#plt.savefig('2018lepto98_gain')
|
|
pylab.show()
|
|
embed()
|
|
|
|
#phase in degree
|
|
# Q10 / conductivity
|
|
# AM-frequency / envelope-frequency scale title? |