Forum: PC-Programmierung
von
Tom (Gast)
09.05.2022 20:42
Hallo zusammen
Ich möchte mittels Python eine einfache Lock-In Verstärkung simulieren.
Dazu generiere ich ein 50Hz Signal welches mit der Zeit abnimmt und ein
100Hz Signal welches mit der Zeit zunimmt, siehe angehängtes Bild der
STFT.
Im 3ten Plot stelle ich den Ausgang der Lock-In Verstärkung dar, die
allerdings keinen Sinn ergibt. Der Ausgang der Lock-In Verstärkung
scheint nicht von der Frequenz "flockin" abzuhängen, wenn ich die
Variable "flockin" (Zeile 28) ändere kommt immer das selbe raus.
Sieht vielleicht jemand den Fehler?
1 from scipy.signal import stft
2 import numpy as np
3 from numpy import sin,pi,exp,sqrt
4 import matplotlib.pyplot as plt
5 from scipy import signal
6
7 dt=1e-3
8 T=10
9 ts=dt*np.arange(0,T/dt)
10 fs=1/dt
11
12 f0=50
13 f1=100
14 tau=4
15 sigs=3*exp(-ts/tau)*sin(2*pi*f0*ts)+0.3*exp(ts/tau)*sin(2*pi*f1*ts)
16
17 n_stft=256
18 fs_,ts_,Cs_=stft(sigs,1/dt,nperseg=n_stft)
19
20 #perform lockin amplification
21 flockin=50
22 s_lockin_d=np.sin(2*pi*ts*flockin)
23 s_lockin_q=np.cos(2*pi*ts*flockin)
24 #compute the measured signal with the d and q component
25 s_prod=np.sqrt((s_lockin_d*sigs)**2+(s_lockin_q*sigs)**2)
26 #filter the signal using a low pass filter
27 flowpass=10
28 flowpass_norm=flowpass/(fs/2)
29 b,a=signal.butter(3,flowpass_norm,'low')
30 lockin_output=signal.filtfilt(b,a,s_prod)
31
32 plt.figure(1)
33 plt.pcolormesh(ts_,fs_,np.abs(Cs_))
34 plt.title("$N_{window}=%d$"%(n_stft))
35 plt.xlabel("Time t [s]")
36 plt.ylabel("DTFT[f](t)")
37
38 plt.figure(2)
39 plt.plot(ts,sigs)
40 plt.xlabel("Time t [s]")
41 plt.ylabel("Signal s(t)")
42
43 plt.figure(3)
44 plt.plot(lockin_output)
von
Tom (Gast)
09.05.2022 21:26
Hmm, das Problem scheint zu sein, dass ich zuerst tiefpassfiltern und
erst naher die Norm der d/q Komponenten berechnen muss, der folgende
Code funktioniert richtig:
1 # -*- coding: utf-8 -*-
2
3 from scipy.signal import stft
4 import numpy as np
5 from numpy import sin,pi,exp,sqrt,square
6 import matplotlib.pyplot as plt
7 from scipy import signal
8
9 dt=1e-3
10 T=10
11 ts=dt*np.arange(0,T/dt)
12 fs=1/dt
13
14 f0=50
15 f1=100
16 tau=4
17 sigs=3*exp(-ts/tau)*sin(2*pi*f0*ts)+0.3*exp(ts/tau)*sin(2*pi*f1*ts)
18
19 n_stft=256
20 fs_,ts_,Cs_=stft(sigs,1/dt,nperseg=n_stft)
21
22 #perform lockin amplification
23 flockin=float(100)
24 s_lockin_d=np.sin(2*pi*ts*flockin)
25 s_lockin_q=np.cos(2*pi*ts*flockin)
26 #compute the measured signal with the d and q component
27 s_prod_d=s_lockin_d*sigs
28 s_prod_q=s_lockin_q*sigs
29 #filter the signal using a low pass filter
30 flowpass=10
31 flowpass_norm=flowpass/(fs/2)
32 b,a=signal.butter(3,flowpass_norm,'low')
33 lockin_output_d=signal.filtfilt(b,a,s_prod_d)
34 lockin_output_q=signal.filtfilt(b,a,s_prod_q)
35
36 plt.figure(1)
37 plt.pcolormesh(ts_,fs_,np.abs(Cs_))
38 plt.title("$N_{window}=%d$"%(n_stft))
39 plt.xlabel("Time t [s]")
40 plt.ylabel("DTFT[f](t)")
41
42 plt.figure(2)
43 plt.plot(ts,sigs)
44 plt.xlabel("Time t [s]")
45 plt.ylabel("Signal s(t)")
46
47 plt.figure(3)
48 plt.plot(ts,np.sqrt(np.square(lockin_output_d))+np.square(lockin_output_q))
49 plt.ylabel("Lockin Output")
50 plt.xlabel("Time t [s]")
Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Thread beobachten
Seitenaufteilung abschalten