2

Following a past question, I'd like to extrapolate a signal like the one below: (red is signal, blue is the the extrapolated )

t1 = 0:0.1:30;
t2 = 0:0.1:40;
s = @(t) sin(0.2*t)+sin(0.5*t+0.05*t.^2); 
figure;
plot(t2,s(t2)); hold on
plot(t1,s(t1));

enter image description here

Linear predictive coding, or AR Burg's method are not successful here. How can you extrapolate that?

Royi
  • 19,608
  • 4
  • 197
  • 238
bla
  • 588
  • 1
  • 4
  • 14
  • 2
    Only commenting as I only have my guess as to what I might do (which I am sure is not optimum). I would low pass filter to estimate and fit the low frequency sinusoidal component, and high pass filter to estimate the high frequency sinusoidal component including it's chirp rate. From this I would then extrapolate the result based on those estimates. – Dan Boschen Sep 13 '18 at 01:19
  • Think of what Burg is doing: it’s an autoregressive method, so it’s assuming that there is a linear phase relationship between samples. If you have an LFM that you’re trying to approximate, what you have is actually quadratic phase information; therefore, Burg isn’t really the right tool here (no auto regressive method is for that matter). Have you looked at any other methods that aren’t AR based? Singular spectrum analysis extrapolation comes to mind – vintagevogue Sep 13 '18 at 01:55
  • yes I'm aware that LPC or Burg's are not appropriate. Can you give more details on the option you mentioned? – bla Sep 13 '18 at 07:52
  • I’ve worked with SSA quite a lot but never for extrapolation, though I’m aware economists use it there (I’ve never had a need). There’s a nice write up on Wikipedia you could/should check out if you’re interested. – vintagevogue Sep 15 '18 at 00:16

1 Answers1

1

Here's the approach I took:

  1. Estimate\find the chirp rate in the data.
  2. stretch (actually compress) the data to "unchirp" it so it becomes periodic.
  3. Use AR Burg.
  4. Undo the stretch based on the chirp rate estimate found in (1)

EDIT: here's some sample code that implements the idea:

dt=0.01;
t = 0:dt:30; % the time vector sampled
t2 = 0:dt:40; % the full time vector to compare the extrap.
s = @(t) sin(0.5*t)+sin(t+0.05*t.^2); % fake data from the question

% 1. find "chirp rate" by finding peaks and fitting a chirp model
[pks locs]=findpeaks(s(t),t);


ModelEqn = 'a*sqrt(x)+b'; % model to fit linear chirp
init_params = [locs(1) 0]'; % for a,b,c of model
f = fit([1:numel(locs)]',locs(:),ModelEqn,'Start', init_params);


% 2.  unchirp (compress) the data to have uniform spaced peaks 
ts=((t-f.b)./f.a).^2;
ti=linspace(ts(1),ts(end),numel(t)); 
s_unchirped=interp1(ts,s(t),ti,'spline');


% 3. extrapolate compressed data via AR Burg
dti=mean(diff(ti));
ti2=ti(1):dti:(ti(1)+2numel(t)dti); % unchirped time to extrap.


N= round(0.75numel(t)); 
a = arburg(s_unchirped, N);
s_extrap = zeros(1, numel(ti2));
s_extrap(1:N) = s_unchirped(1:N); % fill in the known part of the time series
for n=(N+1):numel(ti2)

    s_extrap(n) = -sum(a(2:end) . s_extrap((n-1):-1:(n-N)));
end


% 4. undo compression on s_extrap to transfom back

tiu2=f.a*sqrt(ti2)+f.b; 
s_extrap_unchirped=interp1(tiu2,s_extrap,t2,'spline');


% plot
figure('Position',[100,100,1500,500])
subplot(1,3,1)
plot(t2,s(t2),t(1:N),s(t(1:N)));
legend('Signal ground truth ','Signal to be used for extrap','Location','northoutside')


subplot(1,3,2)
plot(ti, s_unchirped); hold on ; plot(ti2,s_extrap,'r:');
legend('Signal unchirped','Signal unchirped extrap','Location','northoutside')


subplot(1,3,3)
plot(t2,s(t2));hold on; plot(t2,s_extrap_unchirped)
legend('Signal','Signal extrap.','Location','northoutside')

enter image description here

Some thoughts:

  • instead of the linear chirp model used, one can try a more general nonlinear model, for example, a*x.^b+c, and inverse the transform accordingly.

  • AR burg has its limitations here, not sure that is the best way to go with extrapolation of the periodic (unchirped) data, as it needs to sample "enough".

bla
  • 588
  • 1
  • 4
  • 14
  • Would you care to share some more information? The steps are quite clear but there are no results visible. Have you quantified the error (if any), do you have any intuition as to whether this is applicable to other signals too, have you experienced any difficulties and/or caveats while working on it, is there anything else interesting that is worth noting? – ZaellixA Mar 07 '24 at 18:35
  • 1
    sure... I'll edit the answer with some code based on the minimal example in the question... – bla Mar 07 '24 at 19:22
  • @ZaellixA, enjoy the edited answer... – bla Mar 07 '24 at 21:39