Einfache Lock-In Verstärkung in Python - wo liegt der Fehler?

[Tom (Gast)]

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?

from scipy.signal import stft
import numpy as np
from numpy import sin,pi,exp,sqrt
import matplotlib.pyplot as plt
from scipy import signal

dt=1e-3
T=10
ts=dt*np.arange(0,T/dt)
fs=1/dt

f0=50
f1=100
tau=4
sigs=3*exp(-ts/tau)*sin(2*pi*f0*ts)+0.3*exp(ts/tau)*sin(2*pi*f1*ts)

n_stft=256
fs_,ts_,Cs_=stft(sigs,1/dt,nperseg=n_stft)

#perform lockin amplification
flockin=50
s_lockin_d=np.sin(2*pi*ts*flockin)
s_lockin_q=np.cos(2*pi*ts*flockin)
#compute the measured signal with the d and q component
s_prod=np.sqrt((s_lockin_d*sigs)**2+(s_lockin_q*sigs)**2)
#filter the signal using a low pass filter
flowpass=10
flowpass_norm=flowpass/(fs/2)
b,a=signal.butter(3,flowpass_norm,'low')
lockin_output=signal.filtfilt(b,a,s_prod)

plt.figure(1)
plt.pcolormesh(ts_,fs_,np.abs(Cs_))
plt.title("$N_{window}=%d$"%(n_stft))
plt.xlabel("Time t [s]")
plt.ylabel("DTFT[f](t)")

plt.figure(2)
plt.plot(ts,sigs)
plt.xlabel("Time t [s]")
plt.ylabel("Signal s(t)")

plt.figure(3)
plt.plot(lockin_output)

[Tom (Gast)]

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:

# -*- coding: utf-8 -*-

from scipy.signal import stft
import numpy as np
from numpy import sin,pi,exp,sqrt,square
import matplotlib.pyplot as plt
from scipy import signal

dt=1e-3
T=10
ts=dt*np.arange(0,T/dt)
fs=1/dt

f0=50
f1=100
tau=4
sigs=3*exp(-ts/tau)*sin(2*pi*f0*ts)+0.3*exp(ts/tau)*sin(2*pi*f1*ts)

n_stft=256
fs_,ts_,Cs_=stft(sigs,1/dt,nperseg=n_stft)

#perform lockin amplification
flockin=float(100)
s_lockin_d=np.sin(2*pi*ts*flockin)
s_lockin_q=np.cos(2*pi*ts*flockin)
#compute the measured signal with the d and q component
s_prod_d=s_lockin_d*sigs
s_prod_q=s_lockin_q*sigs
#filter the signal using a low pass filter
flowpass=10
flowpass_norm=flowpass/(fs/2)
b,a=signal.butter(3,flowpass_norm,'low')
lockin_output_d=signal.filtfilt(b,a,s_prod_d)
lockin_output_q=signal.filtfilt(b,a,s_prod_q)

plt.figure(1)
plt.pcolormesh(ts_,fs_,np.abs(Cs_))
plt.title("$N_{window}=%d$"%(n_stft))
plt.xlabel("Time t [s]")
plt.ylabel("DTFT[f](t)")

plt.figure(2)
plt.plot(ts,sigs)
plt.xlabel("Time t [s]")
plt.ylabel("Signal s(t)")

plt.figure(3)
plt.plot(ts,np.sqrt(np.square(lockin_output_d))+np.square(lockin_output_q))
plt.ylabel("Lockin Output")
plt.xlabel("Time t [s]")