From 0068d53e7363aeefe83f7ad8e9f012ca8f2b5032 Mon Sep 17 00:00:00 2001
From: Benedikt Seeger <benedikt.seeger@ptb.de>
Date: Wed, 17 May 2023 10:31:39 +0200
Subject: [PATCH] added first version

---
 .gitmodules |  3 +++
 main.py     | 73 ++++++++++++++++++++++++++++++++++++++++++++++-------
 2 files changed, 67 insertions(+), 9 deletions(-)
 create mode 100644 .gitmodules

diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..ebd9f77
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "sinetools"]
+	path = sinetools
+	url = https://gitlab1.ptb.de/TBruns/sinetools.git
diff --git a/main.py b/main.py
index 76d0e8c..81852b6 100644
--- a/main.py
+++ b/main.py
@@ -1,16 +1,71 @@
-# This is a sample Python script.
 
-# Press Umschalt+F10 to execute it or replace it with your code.
-# Press Double Shift to search everywhere for classes, files, tool windows, actions, and settings.
+import numpy as np
+import scipy as sp
+import matplotlib.pyplot as plt
+import sinetools.SineTools as st# This is a sample Python script.
+
+def plot_spectrogram(title, w, fs):
+    ff, tt, Sxx = sp.signal.spectrogram(w, fs=fs, nperseg=64, nfft=128)
+    fig, ax = plt.subplots()
+    ax.pcolormesh(tt, ff[:145], Sxx[:145], cmap='gray_r',
+                  shading='gouraud')
+    ax.set_title(title)
+    fig.show()
+
+def chripphaseFromZero(ts,f0,fstop):
+    #berechnet die thoretirsche phase aktuell gibts noch nen liniearen offset zum sine fit ....
+    f0=f0
+    sweepSeed=(fstop-f0)/ts[-1]
+    phase=2*np.pi*(f0*ts+((sweepSeed*ts*ts)/2))
+    return -1*phase
+def chripLocalSineFitter(t,y,f0,f1,blockSize=8):
+    #fiittet mit der block size sequenziell 3 param sin approxxes und gibt amplitude und phase zurück
+    ft=np.linspace(f0,f1,t.size)
+    abcs=np.zeros((t.size-blockSize,3))
+    fMeanLoacl=np.zeros((t.size-blockSize))
+    for i in range(t.size-blockSize):
+        fMeanLoacl[i]=np.mean(ft[i:i+blockSize])
+        ts=t[i:i+blockSize]
+        ys=y[i:i+blockSize]
+        abcs[i,:]=st.threeparsinefit(ys,ts,fMeanLoacl[i])
+    complexes=abcs[:,0]*1j+abcs[:,1]#komplexe fit ergebnisse ggf. auch nützlich
+    amp= np.abs(complexes)
+    phase=np.angle(complexes)
+    fig,ax=plt.subplots()
+    ax.plot(fMeanLoacl,amp,label="amp")
+    ax.plot(fMeanLoacl,np.unwrap(phase),label='phase')
+    theoreticalChirpPhase=chripphaseFromZero(t,f0,f1)
+    theoreticalChirpPhase=theoreticalChirpPhase[:-blockSize]
+    ax.plot(fMeanLoacl,theoreticalChirpPhase,label="theorie phase")
+    ax.plot(fMeanLoacl,np.unwrap(phase)-theoreticalChirpPhase,label="phase-theorie")
+    ax.legend()
+    fig.show()
+    print("Done")
+    return fMeanLoacl,amp,np.unwrap(phase)
+
 
 
-def print_hi(name):
-    # Use a breakpoint in the code line below to debug your script.
-    print(f'Hi, {name}')  # Press Strg+F8 to toggle the breakpoint.
 
 
-# Press the green button in the gutter to run the script.
 if __name__ == '__main__':
-    print_hi('PyCharm')
+    startTime=0
+    stopTime=0.160
+    points=4096
+    f0=10
+    f1=2500
+    timePoints=np.linspace(startTime,stopTime,num=points)
+    fs=points/(stopTime-startTime)
+
+    chirptimeSignal=sp.signal.chirp(timePoints,f0=f0,f1=f1,t1=timePoints[-1])
+    """
+    fig1,ax1=plt.subplots()
+    plt.title("Linear Chirp time")
+    plt.xlabel('t (sec)')
+    plt.plot(timePoints,chirptimeSignal)
+    plt.show()
+    """
+
+
+    chripLocalSineFitter(timePoints, chirptimeSignal,f0,f1, blockSize=16)
+
 
-# See PyCharm help at https://www.jetbrains.com/help/pycharm/
-- 
GitLab