Method Sinusoidal() imposes an offset

[StefanPofahl]

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))

[Datseris]

Hi,

offset of the decomposed sinusoidal

What does this mean? offset in which respect?

[StefanPofahl]

Offset in vertical direction.
DBG_signal_decompositon.zip
.

[Datseris]

Alright, I guess that's a bug... Feel free to open a PR if you find the source, I don't see myself having time to fix this at this moment :(

[StefanPofahl]

Hi George, I will see, if I will find the time to get into your code. I am not experienced in this kind of action ... Regards, Stefan Am Do., 7. Juli 2022 um 22:15 Uhr schrieb George Datseris < ***@***.***>:

…

Alright, I guess that's a bug... Feel free to open a PR if you find the source, I don't see myself having time to fix this at this moment :( — Reply to this email directly, view it on GitHub <#14 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/APREBE47563LLXKEOMRTUSLVS43EZANCNFSM5242ZAGQ> . You are receiving this because you authored the thread.Message ID: ***@***.***>

-- Stefan Pofahl Zollgasse 5 8020 Graz Österreich Tel.: +43 (316) 33 2001 Mob.: +43(677)63706434