Method Sinusoidal() imposes an offset #14
Description
Hi George,
I see an offset of the decomposed sinusoidal signal, which is in the case of a monochromatic disturbed signal quite easy to cure, if the root cause of this offset in not easy to detect.
If there is more than one sinusoidal signal this trick is not entirely successful :-(
Below my script to reproduce the error.
Regards
Stefan
P.S.:
The plots are interactive, clicking you can show or hide specific lines.
using SignalDecomposition, PlotlyJS, RobustModels
# example script to test the package "SignalDecomposition":
# https://juliadynamics.github.io/SignalDecomposition.jl/dev/
## --- compose signals: clean and disturbed --------------------------------------------------------------------------------
sampling_frequ = 25000 # Sampling frequency
delta_t = 1.0 / sampling_frequ # Time step in sec
n_signal = 5000 # Length of signal
vec_t = collect(range(start=0, step=delta_t, length=n_signal)) # Time vector
t_milli_s = 1000 .* vec_t # time vector mille seconds
frequ_A = 50.0 # [] = Hz
frequ_B = 120.0 # [] = Hz
ampl_A = 0.0 # Amplitude @f1
ampl_B = 1.0 # Amplitude @f2
# Define the frequency domain bin vector / frequency vector:
# It starts with zero and the upper frequency is defined as the so called Nyquist Frequency
# f_Nyquist = 0.5 * Sampling Frequency. vec_f has halve the size of the data points
# in one fft-window.
vec_f = sampling_frequ / n_signal .* collect(range(start=0, step=1, stop=n_signal/2))
# create signal with sinus at frequ_A and frequ_B:
signal_clean = ampl_A .* sin.(2 * pi * frequ_A .* vec_t) + ampl_B .* sin.(2 * pi * frequ_B .* vec_t)
signal_currupted = signal_clean + 2 * randn(size(vec_t))
## --- local functions --------------------------------------------------------------------------
function _PtsXperiods(_frequency::Real, _sampl_rate::Real, _num_periods::Int=1)
if _sampl_rate < 10 * _frequency
error("sampling rate too low, it needs to be at least ten times the investigated frequency!")
end
_pts_of_n_periods = round(Int, _num_periods * _sampl_rate / _frequency)
return _pts_of_n_periods
end
## --- local plot functions ---------------------------------------------------------------------
function plot_time_domaine(_value_vec::Vector{<:Number}, _sampl_frequ::Real, _frequency::Real, _num_periods::Int=1)
_num_pts = _PtsXperiods(_frequency, _sampl_frequ, _num_periods)
_range_ = range(start= 1, step= 1, length= _num_pts)
println("size(_range_): ", size(_range_))
_time_vec = collect(range(start= 0, step= 1 / _sampl_frequ , length= _num_pts))
println("size(_range_): ", size(_range_), ", \t size(_time_vec)", size(_time_vec))
flush(stdout)
line_signal = scatter(; x = _time_vec, y = _value_vec[_range_])
mylayout = Layout(
title_text = "Time Domaine Signal",
xaxis_title_text = "Time / s",
yaxis_title_text = "Time Domaine Signal / -"
)
return plot(line_signal, mylayout)
end
function plot_compare_three_signals(_signal_A::Vector{<:Number}, _signal_B::Vector{<:Number}, _signal_C::Vector{<:Number},
_label_A::AbstractString, _label_B::AbstractString, _label_C::AbstractString, _sampl_frequ::Real, _frequency::Real, _num_periods::Int= 1)
# --- all need to have the same length:
if ~(length(_signal_A) == length(_signal_B) )
error("Missmatch of vector length!")
end
# --- select part of the signal:
_num_pts = _PtsXperiods(_frequency, _sampl_frequ, _num_periods)
# --- build time vector:
_time_vec = collect(range(start= 0, step= 1 / _sampl_frequ , length= _num_pts))
# --- build index vector / index range:
_range = range(start= 1, step= 1, length= _num_pts)
# --- construct plot objects:
line_A = scatter(; x = _time_vec[_range], y = _signal_A[_range], name = _label_A)
line_B = scatter(; x = _time_vec[_range], y = _signal_B[_range], name = _label_B)
line_C = scatter(; x = _time_vec[_range], y = _signal_C[_range], name = _label_C)
_data = [line_A, line_B, line_C]
_layout = Layout(
title_text = String("Three Signals: 1.) " * _label_A * ", 2.) " * _label_B * ", 3.)" * _label_C),
xaxis_title_text = "Time / s",
yaxis_title_text = "Time Domaine Signal / -"
)
return plot(_data, _layout)
end
function plot_filtered_and_input_signal(_signal_clean::Vector{<:Number}, _signal_disturbed::Vector{<:Number}, _signal_filtered::Vector{<:Number}, _signal_residual::Vector{<:Number},
_sampl_frequ::Real, _frequency::Real, _num_periods::Int=1)
# --- all need to have the same length:
if ~(length(_signal_disturbed) == length(_signal_filtered) == length(_signal_residual) )
error("Missmatch of vector length!")
end
# --- select part of the signal:
_num_pts = _PtsXperiods(_frequency, _sampl_frequ, _num_periods)
# --- build time vector:
_time_vec = collect(range(start= 0, step= 1 / _sampl_frequ , length= _num_pts))
# --- build index vector / index range:
_range = range(start= 1, step= 1, length= _num_pts)
# --- construct plot objects:
line_clean = scatter(; x = _time_vec[_range], y = _signal_clean[_range], name = "clean")
line_disturbed = scatter(; x = _time_vec[_range], y = _signal_disturbed[_range], name = "disturbed")
line_filtered = scatter(; x = _time_vec[_range], y = _signal_filtered[_range], name = "filtered")
line_residual = scatter(; x = _time_vec[_range], y = _signal_residual[_range], name = "residual")
_data = [line_clean, line_disturbed, line_filtered, line_residual]
_layout = Layout(
title_text = "Clean, Disturbed, Filtered and Residual Signal",
xaxis_title_text = "Time / s",
yaxis_title_text = "Time Domaine Signal / -"
)
return plot(_data, _layout)
end
## --- main ----------------------------------------------------------------------------------------------------------------
signal_decomposed, signal_residual = decompose(vec_t, signal_currupted, Sinusoidal([frequ_B]))
signal_composed_symetric = signal_decomposed .- mean(signal_decomposed)
# plot_time_domaine(signal_clean, sampling_frequ, frequ_B, 2)
plot_time_domaine(signal_clean, sampling_frequ, frequ_B, 2)
plot_filtered_and_input_signal(signal_clean, signal_currupted, signal_decomposed, signal_residual, sampling_frequ, frequ_B)
plot_compare_three_signals(signal_clean, signal_decomposed, signal_composed_symetric, "original", "decomposed", "symmetric", sampling_frequ, min(frequ_A, frequ_B))