DSP filter problem in the simulation of a Lock-In amplifier

Question

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

Answer 1

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

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

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

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

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

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"}]]}

df["mean noisy mixed"] // plot

lowpass = ExponentialMovingAverage[#, 0.001] &;

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

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

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

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

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

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}