1
$\begingroup$

I wanted To explain to students how a Lock-In amplifier works. The best way is to simulate the behaviour of one. However, in my program, I get unexpected behaviour, which I suspect has to do with numerics.

The way a Lock-In works is that you multiply a sine-wave signal (generator) that is noisy with another sine-wave (mixer) with the same frequency. This gives you a constant at half the amplitude plus one signal at double the frequency (addition of two sin-waves). You then filter this product signal with a low pass filter, which should give you the DC-signal without the noise, which corresponds to the quantity you are interested in. The filter should suppress signals at high frequencies, but it does not in my implementation.

(*frequency*)
υ := 100;
ϕ := 0;
A := 1;
period := 0.1;
step := 0.00001;
noiselevel := 1;
samrate := 100000;
generator[t_] := A*Cos[2*Pi* t *υ];
mixer[t_] := Cos[2 π t (υ + ϕ) ];
(*noisy data*)
dataprop := Table[RandomReal[NormalDistribution[generator[t], 1]], {t, 0, period,
step}]; ListPlot[dataprop]
(*generating a sine-wave table and mixing, eg. multiplying it*)
mixertable := Table[N[mixer[t]], {t, 0, period, step}];
muxi := dataprop*mixertable;
ListPlot[muxi]
(*the filter, a simple low pass one, the smaller a the stronger/slower the filter *)
a := 0.001; 
filtered :=  RecurrenceTable[{smooth[i] == 
 smooth[i - 1] + a*(muxi[[i]] - smooth[i - 1]), smooth[1] == muxi[[1]]},
 smooth, {i, 1, Length[muxi]}];
ListPlot[filtered, PlotRange -> All]

In this plot you can see that despite the filter there is a periodic signal in the output of the filter at 2x the frequency. If you take the Fourier transform you can see the it is at 200Hz.

Periodogram[filtered, FourierParameters -> {1, -1}, 
            SampleRate -> samrate, PlotRange -> {{0, 1000}, {0, 80}}, 
            ScalingFunctions -> "dB"]

If you write the same code in Python the ripple at 200Hz is not there. See here same program in Python

$\endgroup$
4
  • 2
    $\begingroup$ 1) you don't need to use := everywhere, only on the functions. 2) I don't think you need to use RecurrenceTable here. There are errors from muxi[[i]] when I use it - instead use FoldList[#1 + a*(#2 - #1) &, First[muxi], muxi]; $\endgroup$
    – flinty
    Jun 14 '20 at 15:21
  • 2
    $\begingroup$ ^... your smoothing and my Fold are equivalent to ExponentialMovingAverage[muxi, a]. But I'd recommend using LowpassFilter instead. Also it makes sense that the wave would have a spike at 2x the frequency because $\cos^2(x) = (\cos (2 x)+1)/2$. You should account for this effect by writing SampleRate->samrate/2 in your Periodogram. $\endgroup$
    – flinty
    Jun 14 '20 at 16:39
  • $\begingroup$ @flinty, thanks a lot for the comments. The FoldList[] gets rid of the error. I believe that ExponentialMovingAverage[muxi,a] is exactly the same filter. However, the oscillations in the output persist. I believe they should not be there, as in the Python code these oscillations are not there. I ended up trying the Python example, as I did not understand what is going on. $\endgroup$
    – majeriisli
    Jun 14 '20 at 17:39
  • 4
    $\begingroup$ "If you write the same code in Python the ripple at 200Hz is not there. " But it's not the same code at all. The performance of filter is highly influenced by the parameters. Though the frequecy in Python code is also 100, the period there is 1, while you choose period=0.1. Just choose \[Upsilon] := 1000;, the output of Mathematica will be similar to that of Python. BTW, you're generating the noisy data in an unusual way (it doesn't significantly influence the result though), how do you figure out this undocumented syntax of RandomReal (Seems to be equivalent to RandomVariate)? $\endgroup$
    – xzczd
    Jun 15 '20 at 4:21
3
$\begingroup$

As mentioned in the comment above, your code is not the same as the Python code. Just for fun, here's the same (to be precise, do almost the same job) one. I've intentionally mimicked the style of the Python code, on the premise of keeping the code clean.

P.S. The mixing of sine and sin in the Python code is just terrible, so I don't follow that convention.

sinwave[freq_, phase_: 0, Fs_: 10^4] := Sin[2 Pi freq Range[0., Fs - 1]/Fs + phase Degree]
sin = sinwave[100];
meanarray = ConstantArray[Mean@#, Length@#] &;
mean = meanarray@sin;

df = <|"sin" -> sin, "mean" -> mean|>;

plot = ListLinePlot[#, Frame -> True, Axes -> False, 
                       PlotLegends -> Automatic, PlotRange -> All] &;

df[[All, ;; 1000]] // plot

enter image description here

df["sin mixed"] = sin^2;

df["mean mixed"] = df["sin mixed"] // meanarray;

df[[{"sin mixed", "mean mixed"}, ;; 1000]] // plot

enter image description here

df["sin mixed 101"] = sin sinwave[101];
df["mean mixed 101"] = df["sin mixed 101"] // meanarray;
df[[{"sin mixed 101", "mean mixed 101"}]] // plot

enter image description here

noise = RandomReal[{-1, 1}, 10^4];

df["sin noisy"] = 0.1 sin + noise;
df["sin noisy mixed"] = df["sin noisy"] sin;
df["mean noisy mixed"] = df["sin noisy mixed"] // meanarray;
plot /@ {df["sin noisy"], df[[{"sin noisy mixed", "mean noisy mixed"}]]}

enter image description here

df["mean noisy mixed"] // plot

enter image description here

lowpass = ExponentialMovingAverage[#, 0.001] &;

df["sin mixed lp"] = lowpass@df["sin mixed"];
df["sin mixed lp"] // plot

enter image description here

df["sin noisy mixed lp"] = lowpass@df["sin noisy mixed"];
df["sin noisy mixed lp"] // plot

enter image description here

df["sin noisy mixed lp2"] = lowpass@df["sin noisy mixed lp"];
df["sin noisy mixed lp2"] // plot

enter image description here

df["sin phase45 mixed"] = df["sin noisy"] sinwave[100, 45];
df["sin phase90 mixed"] = df["sin noisy"] sinwave[100, 90];

df["sin phase45 mixed lp"] = lowpass@df["sin phase45 mixed"];
df["sin phase90 mixed lp"] = lowpass@df["sin phase90 mixed"];
    
df[[{"sin noisy mixed lp", "sin phase45 mixed lp", "sin phase90 mixed lp"}]] // plot

enter image description here

df[[{"sin noisy mixed lp", "sin phase45 mixed lp", "sin phase90 mixed lp"}, 
    6000 ;;]] // plot

enter image description here

coswave[freq_, phase_: 0, Fs_: 10^4] := Cos[2 Pi freq Range[0., Fs - 1]/Fs + phase Degree]

df["cos noisy mixed"] = df["sin noisy"] coswave[100];
df["cos noisy mixed lp"] = lowpass@df["cos noisy mixed"];
df["noisy quad mag"] = Sqrt[df["cos noisy mixed lp"]^2 + df["sin noisy mixed lp"]^2];
df["noisy quad pha"] = ArcTan[df["cos noisy mixed lp"], df["sin noisy mixed lp"]];
{df["noisy quad mag"][[8000 ;;]] // plot, 
 df["noisy quad pha"][[8000 ;;]] // plot}

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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