Simulation of Lock-In Amplification in Python makes no sense

Question

I would like to simulate a basic lock-in amplification for post-processing some data. I wrote a basic Python script for this purpose but the output of the lock-in amplification does not seem to make sense. I am simulating a 50Hz signal that increases with time and a 100Hz signal that decreases with time, see the following STFT plot:

$s(t)=3*e^{-t/5}\sin(2\pi50t)+0.3e^{t/4}\sin(2\pi100t)$

The output of the lockin-amplification does not make sense and it does not seem to change when I change the lockin-frequency:

Here is my python code:

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(ts,lockin_output)
plt.ylabel("Lockin Output")
plt.xlabel("Time t [s]")

Does anyone know where the issue could be?

Answer 1

Ok, the problems seems to be that I need to first low-pass filter both the d and q components separately before applying the norm operation, here is the corrected code:

# -*- 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]")