06.08

sin_response_fit.py

@@ -1,28 +1,83 @@
 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
 
-data = np.load('files.npy')
+def take_second(elem):
-for d in data:
+    return elem[1]
+
+
+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)
+    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, jar - np.mean(jar))
-    plt.plot(time, y)
+    #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)
 
-    plt.plot(time, sin_response(time, *sinv), label='fit: f=%f, p=%.2f, A=%.2f' % tuple(sinv))
+    phaseshift.append(np.sqrt(sinv[1]**2))
+    gain.append(np.sqrt(sinv[2]**2))
+    amfreq.append(d[1])
+
+
+    #plt.plot(time, sin_response(time, *sinv), label='fit: f=%f, p=%.2f, A=%.2f' % tuple(sinv))
 
     # plt.legend()
-    plt.show()
+    #plt.show()
-embed()
+
+    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)
+
+fig = plt.figure()
+ax = fig.add_subplot(1, 1, 1)
+ax.plot(amf, mgain, 'o')
+ax.set_yscale('log')
+pylab.show()
 
 #betrag von A

sin_response_specto.py

@@ -38,8 +38,13 @@ datasets = [#'2020-06-19-aa',   #-5Hz delta f, horrible fit
             '2020-07-21-al',      #sin
             '2020-07-21-am',
             '2020-07-21-ak',
-            #'2020-07-21-an',
+            '2020-07-21-an',
-            #'2020-07-21-ao'
+            '2020-07-21-ao',
+            '2020-07-22-ai',
+            '2020-07-22-aj',
+            '2020-07-22-ak',
+            '2020-07-22-al',
+            '2020-07-22-am',
             ]
 
 time_all = []
@@ -83,7 +88,9 @@ for idx, dataset in enumerate(datasets):
 
         file_name.append(str(d))
         files.append(file_name)
-
+        if float(amfreq) < 0.01:
+            spec, freqs, times = specgram(dat, Fs=1/dt, detrend='mean', NFFT=nfft, noverlap=nfft * 0.8)
+        else:
             spec, freqs, times = specgram(dat, Fs=1 / dt, detrend='mean', NFFT=nfft, noverlap=nfft * 0.95)
 
         dbspec = 10.0*np.log10(spec)  # in dB