Free-knot spline approximation (BSFK) problem

@BrunoLuong

My data acquisition system produce periodically 1-D measured noised data with the fixed time window length W. I want to produce smoothed data for each window W separately, with specific constraints on continuous (k=2) or smooth (k = 3 or 4) processed signal connections between consecutive time measurement windows. So, for each window W I get finally separate "pp" structure. How to set proper BSFK options setting to fulfil these constraints?

The second question is: Is there any method how to merge separate "pp" structures to one "pp" structure for several processing windows at one?

Add note: May by some processing windows overlap could be required. Do you have any experiance with using BSFK in streaming regime?

[Bruno Luong]

I'm not I fully understand: why can't you fit the whole sets of data by concatenate them together? It wil provide the single pp and not several that you won't bother to merge

[Michal]

Due to the required accuracy of the fit is processing of the whole dataset extremely difficult (time and memory consuming). Typical length of processing window is about 1e4 samples. Signal is sometimes rapidly changing. Nknots ~ 60-100 for window W. Total number of Windows is really very high (~1e3 - 1e4).

I need near real-time processing. So incremental processing scenario is demanded.

[Bruno Luong]

I see.

Can you share few consecutive pp structs (ideally 3-5) where you want to merge?

[Michal]

See attached MAT file (result_4_8.mat) with complete BSFK results for 4 different types of signal and 8 consecutive processing windows.

structure result{SigNum,Wnum} containts all relevant data (pp,x,y, params, options , ...).

SigNum = 1,2,3,4 and Wnum = 1,2, ... 8. Typical elapsed time for 1 processing window is about 35-45secs (8core CPU)

How dangerous is the following warning (for SigNum = 2, Wnum = 5,7)?

Warning: 0 of the 1 requested eigenvalues converged. Eigenvalues that did not converge are NaN. 
> In eigs (line 176)
In BSFK>RestoreSPMat (line 2827)
In BSFK>GaussNewtonStep (line 2306)
In BSFK (line 494)

[Bruno Luong]

This warning is not dangereous in my experience.

[Michal]

Opps ... The attached MAT file is missing!!?? Something is wrong with attachments ...

[Michal]

The MAT file with all relevant data ... 2nd attempt :)

[Bruno Luong]

What I can offer is reduced the data using median filter

data = load('result_4_8.mat')
data = data.result;

[m,n] = size(data);
x = cellfun(@(data) data.x, data, 'unif', 0);
y = cellfun(@(data) data.y, data, 'unif', 0);

j = 1; % select which column
xj = cat(1, x{:,j});
yj = cat(1, y{:,j});

% Assuming the number of data is divisible by 100
reshape_xj = reshape(xj, 100, []);
reshape_yj = reshape(yj, 100, []);
reduced_xj = median(reshape_xj,1);
reduced_yj = median(reshape_yj,1);
pp = BSFK(reduced_xj,reduced_yj, 4);

% Graphic check
xi = linspace(min(xj),max(xj),1000);
yi = ppval(pp, xi);
plot(xj, yj,'c.');
hold on
plot(xi, yi, 'r', 'Linewidth', 2)

[Bruno Luong]

Result with

pp=BSFK(reduced_xj,reduced_yj,2,100);

[Michal]

OK ...

So, you do not see any chance how to connect two separate consecutive splines?

Using shape BSFK constraints, for Example

[Bruno Luong]

BSFK is not supposed to do that. Its a fitting function and it produces a single pp.

You could always connect 2 piecewise functions in C1 manner by a 3rd order polynomial if there is a gap between them, since there are 4 conditions and 4 unknown.

It is likely not produce a nice transition since one pp ignores completely the data next to it.

[Bruno Luong]

You could try to recursively enforce the continuity for function/derivative when you call BSFK on the next interval using the pp of the previous interval, just tell BSFK to have function/derivative of the most left knot (current) = previous pp function/derivative at the right knot (previous).

[Michal]

Yes, something like this I would like to try, but I don't know how exactly implement this type of fix by "shape" and "pntcon" options. This is the main reason why I need your help.

Answer 1

[Bruno Luong]

Reference to BSFK function https://fr.mathworks.com/matlabcentral/fileexchange/25872-free-knot-spline-approximation

Here is the recursive pointwise constraint. You'll see it does the job (zoom in) the transition is not nice

data=load('result_4_8.mat')
data=data.result;
[m,n] = size(data);
x = cellfun(@(data) data.x, data, 'unif', 0);
y = cellfun(@(data) data.y, data, 'unif', 0);
j = 1; % n
close all
figure
hold on
for i = 1:m
    % Normalize data so that dy/dx is comparable to y
    xij = x{i,j}/10000;
    yij = y{i,j}/10000;
    options = struct('lambda', 1e-8);
    if i >= 2
        xleft = pp.breaks(end);
        yleft = ppval(pp,xleft);
        ydleft = ppval(ppder(pp),xleft);
        xij = [xleft; xij];
        yij = [yleft; yij];
        pntcon = struct('p', {0 1}, 'x', {xleft,xleft}, 'v', {yleft ydleft});
        options.pntcon = pntcon;
    end
    pp =BSFK(xij, yij, 4, [], [], options);
    xi = linspace(min(xij),max(xij),1000);
    yi = ppval(pp, xi);
    plot(xij, yij,'c.');
    plot(xi, yi, 'r', 'Linewidth', 2);
    drawnow
end

function ppd = ppder(pp)
ppd = pp;
coefs = ppd.coefs;
n = size(coefs,2);
ppd.coefs = coefs(:,1:n-1).*(n-1:-1:1);
ppd.order = ppd.order-1;
end

[Bruno Luong]

You can play with making data on segment overlapping, etc... based on the syntax I gave you.

Good luck.

[Michal]

Thanks Burno!!!

[Michal]

@Bruno Luong Are there some theoretical numerical complexity and memory requirements estimation for BSFK depending on approximation order (k), number of signal points and number of knots (nknots)?

My current experience with BSFK:

1. CPU time increase at least quadraticaly with number of signal points
2. CPU time increase at least linearly with number of knots (nknots)
3. CPU time increasy at least linearly with approximation polynom order (k)

But all these observation I have made are so far not reliable.

[Bruno Luong]

I haven't not studied the complexity of BSFK.

But it seems to me the linear dependency to number of knots is not quite true, I would say it is more like quadratic.

[Michal]

Of course, this is my mistake ... you are right! CPU time dependency to number of knots is ~ quadratic!!!