# Detection of constant, weak signals
**URL:** https://discourse.julialang.org/t/detection-of-constant-weak-signals/101022
**Category:** Signal and Image Processing
**Tags:** question
**Created:** [June 30, 2023, 3:22pm UTC](https://discourse.julialang.org/t/detection-of-constant-weak-signals/101022 "2023-06-30T15:22:23Z")
**Posts on this page:** 12
**Page:** 1
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [June 30, 2023, 3:22pm UTC](https://discourse.julialang.org/t/detection-of-constant-weak-signals/101022/1 "2023-06-30T15:22:23Z")
I have the following problem:
I have signal u(t), which can be sinus function, but I could
also use any other function as exitation signal.
I have a system which has a gain k between -1 and 1 and adds a
strong noise with a known spectrum (not white noise).
I measure the resulting signal v(t)=k \* u(t) + noise(t).
The noise has an amplitude which is about 100 times of the amplitude
of u(t). I want to determine k as fast and accurate as possible.
How can I do this?
---
### Author:  [@mbaz](https://discourse.julialang.org/u/mbaz)
#### Post date: [June 30, 2023, 3:40pm UTC](https://discourse.julialang.org/t/detection-of-constant-weak-signals/101022/2 "2023-06-30T15:40:05Z")
I can think of a couple of different approaches. One would be to to apply a whitening filter (to remove the correlation in the noise samples) with unity gain at the frequency of `u(t)`, followed by a bandpass filter and then correlation between the whitened signal and `u(t)`. It should be possible to infer `k` from the correlation.
Another approach is to take advange of the fact that you know the frequency and phase of `u(t)` to calculate a 3-parameter LS sinusoidal regression on v(t) (maybe after whitening).
Finally, since the noise is not white, you may want to choose the frequency of `u(t)` to be in the “less noisy” band of the system. Apply `u(t)`, then a very narrow bandpass filter. Again, you may be able to infer `k` from the magnitude of the DTFT of the output.
In any case you’re likely to need a very long measurment time, to let `u(t)` have enough energy to be detected.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [June 30, 2023, 3:43pm UTC](https://discourse.julialang.org/t/detection-of-constant-weak-signals/101022/3 "2023-06-30T15:43:15Z")
> [@ufechner7](#):
>
> a strong noise with a known spectrum […] I want to determine k as fast and accurate as possible.
In the Fourier domain, this is white noise with a known diagonal scaling. In this case isn’t diagonally weighted least-squares (in the Fourier domain) statistically optimal (for linear estimators) by the [Gauss–Markov theorem](https://en.wikipedia.org/wiki/Gauss%E2%80%93Markov_theorem)?
Since there is only one unknown, least-squares is simply a dot product (a 1d projection): FFT your data v(t) (maybe with a window function to ameliorate edge effects) and your function u(t) to obtain vectors \hat{v} and \hat{u}, and then the weighted least-squares k is simply `k = dot(û, S, v̂) / dot(û, S, û)` with `S = Diagonal(inv.(noisespectrum))`.
(In the time domain, this is essentially equivalent to a correlation coefficient with a whitening filter as suggested by @mbaz.)
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [June 30, 2023, 3:47pm UTC](https://discourse.julialang.org/t/detection-of-constant-weak-signals/101022/4 "2023-06-30T15:47:54Z")
What I do in the moment is to use lock-in detection, see [Principles of Lock-in Detection | Zurich Instruments](https://www.zhinst.com/europe/en/resources/principles-of-lock-in-detection) . It is somehow working, but I am looking for a better optimized approach. Lock-in detection does NOT assume a constant signal, and you need to define a border frequency for the low pass filter, which is not what I want. What I want is a result that gives me a. an estimate of k and b. something like a variance or other error property.
One hour of measurement is fine, longer might also be acceptable.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [June 30, 2023, 3:53pm UTC](https://discourse.julialang.org/t/detection-of-constant-weak-signals/101022/5 "2023-06-30T15:53:32Z")
> [@ufechner7](#):
>
> something like a variance or other error property.
Least-squares estimation also has a simple formula for the variance of k, especially for weighted least squares with \operatorname{Var}(\hat{v}) = E[\hat{v}\hat{v}^\*] - E[\hat{v}] E[\hat{v}]^\* = E[\hat{n}\hat{n}^\*] = S^{-1} as above. In this case I believe it works out to simply \operatorname{Var}(k) = \frac{1}{\hat{u}^\*S \hat{u}}, i.e. `1/dot(û, S, û)`. (But apply this formula you need to pay careful attention to the scaling of your FFT, to ensure that `inv(S)` is really the covariance matrix of `v̂` and not some multiple thereof. To be safe, especially in case you got your noise spectrum S slightly wrong, I would use something like `length(v̂) * var(v̂ - k*û) / (tr(inv(S)) * dot(û, S, û))` as your estimate for the variance of k.)
---
### Author:  [@mgkuhn](https://discourse.julialang.org/u/mgkuhn)
#### Post date: [July 1, 2023, 3:39pm UTC](https://discourse.julialang.org/t/detection-of-constant-weak-signals/101022/6 "2023-07-01T15:39:14Z")
Assuming your noise is zero mean and wide-sense stationary, meaning that all you can possibly know about the noise is described by its power spectrum (or equivalently its auto-correlation function), I’d calculate the estimate \hat k = \frac{\int U^\*(f)\cdot V(f)\text{d}f}{\int U^2(f)\text{d}f} while choosing u(t) with a view to minimizing \hat n = \frac{\int U^\*(f)\cdot \text{Noise}(f)\text{d}f}{\int\_0^T U^2(f)\text{d}f}, for example by constructing u(t) such that its power spectrum |U(f)|^2 is proportional to the inverse of the power spectrum of the noise. Or you could make u(t) a sine wave with a frequency chosen at the minimum of the power spectrum of the noise. (The latter is more or less just the idea behind a lock-in amplifier.)
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [July 1, 2023, 5:41pm UTC](https://discourse.julialang.org/t/detection-of-constant-weak-signals/101022/7 "2023-07-01T17:41:27Z")
> [@mgkuhn](#):
>
> I’d calculate the estimate …
The formula you propose is ordinary least squares, whereas the optimal linear estimator should be least squares weighted inversely by the noise spectrum as noted above. Also, I don’t think the OP can choose u(t)?
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [July 1, 2023, 10:44pm UTC](https://discourse.julialang.org/t/detection-of-constant-weak-signals/101022/8 "2023-07-01T22:44:55Z")
> [@ufechner7](#):
>
> strong noise with a known spectrum
Please provide an example.
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [July 2, 2023, 9:54am UTC](https://discourse.julialang.org/t/detection-of-constant-weak-signals/101022/9 "2023-07-02T09:54:12Z")
Example spectrum, exitation at 0.02 Hz.
I can see already that it would be better to choose a different exitation frequency.
I am trying to write a MWE so that you can have a look at my code, but that is not
so easy…
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [July 2, 2023, 11:13am UTC](https://discourse.julialang.org/t/detection-of-constant-weak-signals/101022/10 "2023-07-02T11:13:28Z")
OK, here my MWE:
```julia
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
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?
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [July 3, 2023, 11:03am UTC](https://discourse.julialang.org/t/detection-of-constant-weak-signals/101022/11 "2023-07-03T11:03:18Z")
Your approach seems good to me. Based on the information provided, I can’t think of anything better.
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [July 3, 2023, 11:12am UTC](https://discourse.julialang.org/t/detection-of-constant-weak-signals/101022/12 "2023-07-03T11:12:50Z")
Well, in literature you can find at least 6 methods for weak signal detection: [https://www.atlantis-press.com/article/25893673.pdf](https://www.atlantis-press.com/article/25893673.pdf)
FFT is one, but others are Wavelet Analysis, adaptive filtering and the Stochastic Resonance Method.
I will try a wavelet analysis next…
As clarification: My noise is just the derivative of the wind speed, to be precise the rotor effective wind speed. The noise spectrum I provided here is very much simplified and depends on the average wind speed, the location and the weather…