OK, here my MWE:

```
using PythonPlot, FFTW, MAT, Statistics
N = 100000 # number of samples
dt = 0.5 # sampling time
ex::Float64 = -0.02 # amplitude of exitation signal, 2% of the noise amplitude
f_ex::Float64 = 0.02 # exitation frequency in Hz
const WINDOW = 20 # number of frequencies to average for the error estimate
LOGPLOT = true
mutable struct Spectrum
spec
freqs
end
function low_pass(signal, gain, initial=0.0)
res = similar(signal)
last_in = signal[begin]
last_out = initial
i = 1
for value in signal
res[i] = last_out + (value-last_out)*gain*dt
last_in = value
last_out = res[i]
i += 1
end
res
end
# https://en.wikipedia.org/wiki/High-pass_filter
function high_pass(signal, alpha)
res = similar(signal)
last_in = signal[begin]
last_out = 0.0
i = 1
for value in signal
res[i] = alpha * last_out + alpha * (value-last_in)
last_in = value
last_out = res[i]
i += 1
end
res
end
function differentiate(signal)
res = similar(signal)
last_value = signal[begin]
i = 1
for value in signal
res[i] = value-last_value
last_value = value
i += 1
end
res
end
function real_part!(spec::Spectrum)
len = length(spec.freqs)
first = 3+div(len, 2)
spec.spec = spec.spec[first:end]
spec.freqs = spec.freqs[first:end]
end
function plot_fft(signal)
global spec, N
signal .-= mean(signal)
N = length(signal)
tmax = (N-1) * dt
t = 0:dt:tmax
# Fourier Transform of it
F = fft(signal) |> fftshift
freqs = fftfreq(length(t), 1.0/dt) |> fftshift
spec = Spectrum(abs.(F), freqs)
spec_phase = Spectrum(angle.(F), freqs)
real_part!(spec)
real_part!(spec_phase)
index = findfirst(==(f_ex), spec.freqs)
k_est = -spec.spec[index] * sign(spec_phase.spec[index])/N*2
# calculate error estimate
err_est = 0.0
for i in -WINDOWĂ·2:WINDOWĂ·2
if i != 0
err_est += spec.spec[index+i]
end
end
err_est /= WINDOW
k_est_error = err_est / N * 2
println("amplitude: $(spec.spec[index]), err_est: $(err_est)")
println("phase: $(spec_phase.spec[index])")
println("k_est: $(round(k_est, digits=3)), k_est_error: $(round(k_est_error, digits=3)) ")
if LOGPLOT
loglog(spec.freqs, spec.spec, label = "Spectrum")
else
plot(spec.freqs, spec.spec, label = "Spectrum")
xlim(f_ex/2, f_ex*1.5)
ylim(0, 1500)
end
grid(true)
nothing
end
# create white noise of the desired length
rews = (rand(N).-0.5)*10.0
# apply some filters to achieve the desired spectrum
rews_filt1 = low_pass(rews, 0.03)
rews_filt = low_pass(rews_filt1, 0.2)
rews_filt = rews_filt .+= 4*high_pass(rews_filt, 0.4)
rews_filt = rews_filt .+= 2*high_pass(rews_filt, 0.1)
tmax = (N-1) * dt
t = 0:dt:tmax
ex_in = ex * sin.(2Ď€ * t * f_ex)
diff_rews = differentiate(rews_filt) ./ dt # sampling frequency is 2 Hz
amp_old = (maximum(diff_rews)-minimum(diff_rews)) / 2.0
diff_rews ./= amp_old # make sure the amplitude of the noise is one
sum_in = ex_in .+ diff_rews
plot_fft(sum_in)
title("Measured signal")
nothing
```

Sorry for the lengthâ€¦

Output:

```
julia> include("src/mwe.jl")
amplitude: 1119.2794208164655, err_est: 224.94989050794956
phase: 1.5655435549107437
k_est: -0.022, k_est_error: 0.004
```

If you change the value for ex, then k_est should change accordinglyâ€¦ If you run the program multiple times you will get different estimates.

Zoom in around the exitation frequency:

I calculate an error estimate by averaging the amplitude for 10 frequencies above and below the exitation signalâ€¦

Is that a valid approach?