0
$\begingroup$

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)$

enter image description here

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

enter image description here

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?

$\endgroup$

1 Answer 1

1
$\begingroup$

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]")
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.