From 183c95cb8b05a302d9e915363452c0d9ec8adeb5 Mon Sep 17 00:00:00 2001
From: xaver <xaver.roos@t-online.de>
Date: Wed, 5 Aug 2020 17:16:49 +0200
Subject: [PATCH] 05.08

---
 scratch.py             |  4 +--
 sin_response_fit.py    | 28 ++++++++++++++++++++
 sin_response_specto.py | 59 +++++++++++++++++++-----------------------
 3 files changed, 57 insertions(+), 34 deletions(-)

diff --git a/scratch.py b/scratch.py
index 1e715c2..bd33e4a 100644
--- a/scratch.py
+++ b/scratch.py
@@ -44,5 +44,5 @@ for d in range(int(len(data)/10)):
 g = [1.2917623576698833, -5.479055166593157, -2.689492238578325, -0.11604244418416806, -0.05353823781665627]
 a = [0.2, 0.002, 0.02, 0.5, 1.0]
 
-plt.plot(a, g)
-plt.show()
+np.save('g.npy', g)
+print(np.load('g.npy'))
diff --git a/sin_response_fit.py b/sin_response_fit.py
index e69de29..ea2122e 100644
--- a/sin_response_fit.py
+++ b/sin_response_fit.py
@@ -0,0 +1,28 @@
+from scipy import signal
+import matplotlib.pyplot as plt
+import numpy as np
+from IPython import embed
+from scipy.optimize import curve_fit
+from jar_functions import sin_response
+
+data = np.load('files.npy')
+for d in data:
+    dd = list(d)
+    jar = np.load('%s.npy' %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, 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))
+
+    # plt.legend()
+    plt.show()
+embed()
+
+#betrag von A
\ No newline at end of file
diff --git a/sin_response_specto.py b/sin_response_specto.py
index 8e40c87..49203ba 100644
--- a/sin_response_specto.py
+++ b/sin_response_specto.py
@@ -47,19 +47,13 @@ freq_all = []
 
 amfrequencies = []
 gains = []
+files = []
+
 
 ID = []
 col = ['dimgrey', 'grey', 'darkgrey', 'silver', 'lightgrey', 'gainsboro', 'whitesmoke']
 labels = zip(ID, datasets)
 
-for infodataset in datasets:
-    infodataset = os.path.join(base_path, infodataset, 'info.dat')
-
-    i = parse_infodataset(infodataset)
-    identifier = i[0]
-    if not identifier[1:-2] in ID:
-        ID.append(identifier[1:-2])
-
 
 for idx, dataset in enumerate(datasets):
     datapath = os.path.join(base_path, dataset, '%s.nix' % dataset)
@@ -71,16 +65,26 @@ for idx, dataset in enumerate(datasets):
     nfft = 2**17
 
     for d, dat in enumerate(data):
-        amfreq = float(import_amfreq(datapath))
+        file_name = []
+
+
+        for infodataset in datasets:
+            infodataset = os.path.join(base_path, infodataset, 'info.dat')
+
+            i = parse_infodataset(infodataset)
+            identifier = i[0]
+            if not identifier[1:-2] in ID:
+                ID.append(identifier[1:-2])
+        file_name.append(ID[0])
+
+        amfreq = import_amfreq(datapath)
         print(amfreq)
-        amfrequencies.append(amfreq)
+        file_name.append(str(amfreq))
 
-        file_name = []
-        file_name.append(ID)
-        file_name.append(amfreq)
-        file_name.append(d)
+        file_name.append(str(d))
+        files.append(file_name)
 
-        spec, freqs, times = specgram(dat, Fs=1/dt, detrend='mean', NFFT=nfft, noverlap=nfft*0.6)
+        spec, freqs, times = specgram(dat, Fs=1/dt, detrend='mean', NFFT=nfft, noverlap=nfft*0.95)
 
         dbspec = 10.0*np.log10(spec)  # in dB
         power = dbspec[:, 50]
@@ -106,30 +110,21 @@ for idx, dataset in enumerate(datasets):
         jm = jar4 - np.mean(jar4)
 
         cut_times = times[:len(jar4)]
-        #np.save('%s.npy' % file_name, jar4)
+        np.save('time: %s.npy' % file_name, cut_times)
+        np.save('%s.npy' % file_name, jar4)
 
         plt.plot(cut_times, jm, '-k')
 
-        cf, ct = mean_noise_cut(jar4, cut_times, n = int(round(len(jar4)/((times[-1] - times [0]) * amfreq))))
-
+        #cf, ct = mean_noise_cut(jar4, cut_times, n = int(round(len(jar4)/((times[-1] - times [0]) * amfreq))))
         #plt.plot(ct, cf, '-k')
 
         #plt.imshow(spec4, cmap='jet', origin='lower', extent=(times[0], times[-1], lim0, lim1), aspect='auto', vmin=-80, vmax=-10)
-        #np.save( , spec4)
-
-        embed()
-        b, a = signal.butter(4, 0.01 / 10000, 'high', analog=True)
-        y = signal.filtfilt(b, a, jm)
-        sinv, sinc = curve_fit(sin_response, cut_times, jm, [amfreq, 2, 0.5])
-        print('frequency, phaseshift, amplitude:', sinv)
-        gains.append(sinv[2])
-
-        plt.plot(cut_times, sin_response(cut_times, *sinv), label='fit: f=%f, p=%.2f, A=%.2f' % tuple(sinv))
 
         #plt.legend()
         #plt.ylim(lim0, lim1)
-        plt.legend()
-        plt.show()
+        #plt.legend()
+        #plt.show()
+        np.save('files.npy', files)
+#embed()
 
-        #print(np.load('%s.npy' % file_name))
-embed()
\ No newline at end of file
+# running average over on AM-period?
\ No newline at end of file