interp1 not working as expected and returning the same negative index for multiple query points (sample code and sample data included)

Hi!

I'm using interp1 to estimate indices that correspond to a magnitude threshold before and after a peak I've identified on a curve. My approach was to segment the curve into a pre-peak dataset, and a post-peak data set and use interp1 on each dataset to find where the curve crosses this threshold on either side of the peak. I've included a sample figure to help visualize this. I've also included some sample code and sample data in case that helps.

So, using interp1, I would like to estimate the exact index at which the data passes this magnitude threshold of -20.

Here is the code I've implemented:

target_val = -20; %some magnitude threshold value

[peak_val, peak_idx] = max(data);

%segmenting the data into pre- and post-peak segments
prePeak_data = data(1:peak_idx);
postPeak_data = data(peak_idx+1:end);

[interp_x, interp_index] = unique(prePeak_data); %using unique() here because I got the "grid vectors must contain unique points" error
prePeak_target_idx = interp1(interp_x, prePeak_data(interp_index), target_val);
[interp_x, interp_index] = unique(postPeak_data);
postPeak_target_idx = interp1(interp_x, postPeak_data(interp_index), target_val);

%also tried re-ordering the variables in interp1() as follows [I think this is the correct implementation but also doesn't work]
[interp_x, interp_index] = unique(prePeak_data); %using unique() here because I got the "grid vectors must contain unique points" error
prePeak_target_idx = interp1(prePeak_data(interp_index), interp_x, target_val);
[interp_x, interp_index] = unique(postPeak_data);
postPeak_target_idx = interp1(postPeak_data(interp_index), interp_x, target_val);

%also tried implementing without using unique()
prePeak_data_x = 1:1:length(prePeak_data);
postPeak_data_x = 1:1:length(postPeak_data);
prePeak_target_idx = interp1(prePeak_data, prePeak_data_x, target_val);
postPeak_target_idx = interp1(postPeak_data, postPeak_data_x, target_val);

Here's the faulty outcome where the indices seem to overlap at -20 (which is also the target magnitude value, so assuming this is not a coincidence). Can confirm that the outcome changes to match whatever target_val I use.

I feel really dumb, thank you in advance for your help! I'm particularly interested in how to get this interp1 to work, but I'll also welcome other suggestions on how to obtain these indices.

longtime dweller, but rare poster. Would defo appreciate constructive feedback on posting norms if applicable

Answer 1

[Star Strider]

Another approach —

DL = load('mathworks help workspace_interp1_20210426.mat');
data = DL.data;

L = numel(data);
samples = linspace(1, L, L);

target_val = -20; %some magnitude threshold value

[peak_val, peak_idx] = max(data);

%segmenting the data into pre- and post-peak segments
prePeak_idx = (1:peak_idx);
postPeak_idx = (peak_idx+1:L);

prePeak_int = find(diff(sign(data(prePeak_idx)-target_val)));                       % Approximate Intercept Index
postPeak_int = find(diff(sign(data(postPeak_idx)-target_val)));                     % Approximate Intercept Index

p1 = polyfit(data(prePeak_int(1)+(-1:1)), samples(prePeak_int(1)+(-1:1)), 1); 
prePeak_x = polyval(p1,target_val);
p2 = polyfit(data(postPeak_int(1)+(-1:1)), samples(postPeak_int(1)+(-1:1)), 1); 
postPeak_x = polyval(p2,target_val);

figure
semilogx(samples(prePeak_idx), data(prePeak_idx),'-b')
hold on
plot(samples(postPeak_idx), data(postPeak_idx),'-b')
plot(prePeak_x, target_val,'xr')
plot(postPeak_x, target_val,'xr')
hold off
grid
xlabel('Samples')
ylabel('Magnitude')
text(prePeak_x, target_val, sprintf('\\leftarrow x = %.2f',prePeak_x), 'Horiz','left','Vert','middle')
text(postPeak_x, target_val, sprintf('x = %.2f \\rightarrow ',postPeak_x), 'Horiz','right','Vert','middle')

producing —

With even dlightly noisy data, the interp1 function will not work correctly, so I opted for linear regression instead.

I cannot load the .mat file in the onlilne application so I cannot run it here.

.

[Karam Elabd]

Thank you so much!! This seems to be doing an ok job of estimating the points of interest. When I ran the code across the remainder of my dataset, it appeared like the index approximations deviated a fair bit away from the curve. see below for all samples in a subset and a close-up view showing the misalignment.

I wrote it off to the data being too noisy but didn't investigate much further. I resorted to marinating a little longer on the interp1 function and I found my problem. It was silly as expected - I believe I didn't pass the interp1 function the correct y-variables, and so my function was returning NaNs because the query point was out of range (classic, right?). Here's the product after some progress:

I definitely implemented some of your coding style and the semilogx respresentation is going to be neater for presenting this data!

Thanks again Star Strider and Mathieu

[Star Strider]

As always, my pleasure!

The data are much noisier than I initially realised.

That problem can be solved by increassing the range of indices on the descending part of the data so the noise effect is significantly reduced —

p2 = polyfit(data(postPeak_int(1)+(-9:9)), samples(postPeak_int(1)+(-9:9)), 1); 

Zooming in on the plot of the provided data, after having estimated the data using this expanded version of ‘p2’ and plotting the data points as well as the lines connecting them (with the '.-b' line style) —

So there are actually three potential crossing-points at -20. The linear regression approach has selected the most likely one here, although it could be anywhere between the first and last crossings and still be in the acceptable range.

Note — Because of the extreme slope and many fewer data in the ‘p1’ calculation, it should be kept as it is, and not changed.

.

[Karam Elabd]

Ahh! that totally makes sense. I should've thought to try that. Yea the data is rather noisy, but I've been instructed not to smooth it out. Thanks again Star Strider!

[Star Strider]

As always, my pleasure!

Answer 2

[Mathieu NOE]

hello

why do not use this crossing function to get the positive slope and negative slope crossing point to a given value ?

load('mathworks help workspace_interp1_20210426.mat');
t = 1:length(data);

threshold = -20; % 

[t0_pos,s0_pos,t0_neg,s0_neg]= crossing_V7(data,t,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 
% ind => time index (samples)
% t0 => corresponding time (x) values 
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"

figure(1)
plot(t,data,t0_pos,s0_pos,'+r',t0_neg,s0_neg,'+g','linewidth',2,'markersize',12);grid on
legend('data','positive slope crossing points','negative slope crossing points');

[Karam Elabd]

Thanks so much Mathieu! I've added this function to my arsenal :)

Cheers!

[Mathieu NOE]

you're welcome !