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?
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.
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?
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.)
What I do in the moment is to use lock-in detection, see Principles of Lock-in Detection | Zurich Instruments . 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.
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.)
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.)
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)?
Please provide an example.
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…
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?
Your approach seems good to me. Based on the information provided, I can’t think of anything better.
Well, in literature you can find at least 6 methods for weak signal detection: 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…