Using LPV spectral decompose

[Datseris]

Better than SinuisoidalFit, see chapter 10 of this: https://lup.lub.lu.se/search/ws/files/55399444/phdthesis.pdf

Julia package already exists: https://github.com/baggepinnen/LPVSpectral.jl

[Datseris]

@baggepinnen on Slack you mentioned how this method is superior to a sinusoidal fit and you cited your thesis and this package. However, it is not clear to me how I can use this package to replace a sinuisoidal fit. It is also not clear to me how it fits in signal separation.

I have to use ls_spectral with my signal s and time axis t. A crucial point is that t does not have to be equispaced, and you said on Slack that this is fine. Here is a simple example:

t = cumsum(rand(1000)/2) # non-equispaced timeaxis
periodic = @. 4 + 7.2cos(t) + 5.6cos(2t + 3π/5)
 x, f = ls_spectral(t, periodic .- mean(periodic))
plot(f, abs2.(x))

this doesn't look at all like the spectrum of 2 cosines. What am I doing wrong?

[baggepinnen]

t and y are flipped, try ls_spectral(y,t,f) instead, and selecet the frequencies f to something reasonable for your application. You may also try a windowed version

x, f = ls_windowpsd(periodic, t, 0:0.01:(2pi))
plot(f, abs2.(x))

[Datseris]

Thanks, that was pretty stupid of me... After changing the order and calling

x, f = ls_windowpsd(periodic, t, 0:0.01:(2pi))
plot(f, abs2.(x))

I get something much more reasonable:

However, I must admit I need still help to conceptually understand how I can add this to SignalDecomposition.jl.

1. The peaks appear at frequencies 0.15 and 0.3. I can't seem to make sense of these numbers given my time axis and/or the way I created periodic...
2. What I want to create is the following method: find cosine functions using your method, and then remove them from the original signal. The first step is finding the fourier components of the given frequencies. This is a simple findfirst at the f array.
3. Second step is translating this information to an amplitude and phase of found cosine. I admit I am clueless on how to do this.
4. After I've found the phase and amplitude A, φ, I simply subtract @. A * cos(ω*t + φ) from my signal, with ω the user specified frequency.

I am very unfamiliar with these aspects of signal processing, so I admit it is very confusing for me how to associate the user provided frequency, which in this example is given in inverse units of the actual time axis t, and the resulting frequencies f of ls_windopsd.

If you can provide any help that would be great, because I can incorporate this method into SignalDecomposition.jl before announcing it on discourse.

[baggepinnen]

# 1
julia> 1/(2pi)
0.1591549430918953

# 3
amplitude(x) = abs(x)
phase(x) = angle(x)

You can look at how the functoin get_fourier_regressor is implemented

[Datseris]

Thanks for your support. I'm probably still doing something wrong, because your suggestion doesn't give the correct results:

using LPVSpectral, Statistics, PyPlot, Random

Random.seed!(124541)

tu = cumsum(rand(1000)/2) # non-equispaced timeaxis
pu = @. 4 + 7.2cos(tu) + 5.6cos(2tu + 3π/5)

fs = [1/2π, 1/π]
μ = mean(pu)
x, f = ls_spectral(tu, pu .- μ, fs)

As = abs.(x)
φs = angle.(x)

periodic_component = zero(pu)
for i in 1:length(fs)
    @. periodic_component += As[i]*cos(2π*fs[i]*tu + φs[i])
end

plot(tu, pu)
plot(tu, periodic_component .+ μ)

Is it incorrect to only use this method at only some given frequencies?

[baggepinnen]

Once again, you have flipped the time and signal arguments

[Datseris]

God damn it. Sorry!

The amplitudes are a factor of 2 off when I do it correctly:

[Datseris]

Indeed, doing abs(x)./2 gives a perfect match. Weird.

[baggepinnen]

That factor is probably there to match the output of spectral estimates from DSP. The definition of the amplitude varies between different implementations and textbooks, and they all put parts of 2pi on different places.

[baggepinnen]

https://github.com/baggepinnen/LPVSpectral.jl/blob/c19ed6e482463771a8e7adb59c53b2a940503102/src/lsfft.jl#L35

[Datseris]

Thanks so much, I've got everything I need now. I only have a final question, not so much about coding but about the science behind the method: when doing standard Fourier methods, people advised me to subtract the mean of the signal before applying the methods, because this leads to more accurate results. is this also the case for ls_spectral?

[baggepinnen]

If 0 is in the frequency vector the mean will be automatically estimated, if not, it's best to remove it.

[Datseris]

Hi @baggepinnen , I am faced with a conceptual difficulty when considering whether to remove the mean or not:

I have input signal s. If the seasonal component x of my signal doesn't fully "close" (e.g. I have 1.5 periods), then the seasonal part also has a nonzero mean value. When I am subtracting the mean of my signal, in essence I am also subtracting this mean value.

Does it not matter at all if I do this?

And then, I used to simply re-add the original mean value mean(s) into x, and take their difference as the residual. But of course, now already x has a non-zero mean value.

[Datseris]

Perhaps the best way is to always include 0 frequency, and let the algorithm take care of it?

[baggepinnen]

You are correct, removing the mean in advance is only asymptotically equivalent, including the zero frequency will always do the right thing.