loop error with "unique" and "interp1"

Dear Matlab-Experts,

trying to find the 50% threshold for a sigmoid shaped curve using interp1 for a set of subjects. Managed to do single subjects here:

load cgft
load cgsx
size(cgft)

ans = 1×2

53 101

size(cgsx)

ans = 1×2

8 101

nsub = 8;
[~, ind1] = unique(cgft); % ind = index of first occurrence of a repeated value

ind1 = 5215×1

16 69 122 175 228 281 334 387 440 493

yc(nsub, :) = interp1(cgft(ind1), cgsx(ind1), .50, 'linear','extrap');

Index exceeds the number of array elements. Index must not exceed 808.

But my loop for all 53 subject fails:

% for all subjects:
load cgft
load cgsx
subjects = [1:53];
nbsubjects = length(subjects);

for nsubs = 1:nbsubjects
    nsub = subjects(nsubs);
    [~, ind1] = unique(cgft(nsub, :));
    yc(nsub, :) = interp1(cgft(ind1(nsub, :)), cgsx(ind1(nsub, :)), .50, 'linear','extrap');
end

Would anyone know why the loop doesn't run please? (I'm a matlab dabbler...not a regular user...)

(Error using matlab.internal.math.interp1

Interpolation requires at least two sample points for each grid dimension.)

[Phillip]

upadted cgft file here.

[Torsten]

Both arrays must be of the same size, but they aren't (see above).

[Phillip]

Hi Torsten, thank you for the response. Even when cgsx and cgft are the same size, it doesn't run (only subject by subject but not the loop). The reason why I have to go with "unique" is because of rogue subject 8 that has a very steep S-curve (in red, ignore blue - Capture.jpg).

Answer 1

[Star Strider]

To avoid all the problems about having unique values in one of the vectors, select a narrow range for the interpolation. That (nearly always, in most instances, and always here) prevents that problem.

Try this —

LD1 = load('cgsx.mat');
cgxs = LD1.cgsx;
LD2 = load('cgft.mat');
cgft = LD2.cgft;

for k = 1:size(cgxs,1)
    % qv = cgft(k,:)
    % Q = [max(cgft(k,:)); min(cgft(k,:))]
    val50(k) = 0.5*(max(cgft(k,:))-min(cgft(k,:)))+min(cgft(k,:));              % Detect 50% Of The Range Of A Specific 'cgft' Vector
    ix50 = find(diff(sign(cgft(k,:)-val50(k))));                                % Approximate Index Of That Value
    idxrng = max(1,ix50-1) : min(numel(cgxs(k,:)),ix50+1);                      % Index Range For Interpolation
    cgxs50(k) = interp1(cgft(k,idxrng), cgxs(k,idxrng), val50(k));              % Interpolate On That Range
end

midrange_and_cgxs_values = [val50; cgxs50]

midrange_and_cgxs_values = 2×53

0.5576 0.5230 0.6505 0.4590 0.5960 0.5204 0.5932 0.6040 0.5709 0.5737 0.5642 0.6004 0.4945 0.4715 0.5386 0.4730 0.5525 0.4805 0.5025 0.5074 0.5127 0.6481 0.5234 0.6650 0.4954 0.6038 0.5425 0.6186 0.5018 0.5007 -4.6388 -5.5379 -5.4702 -3.4883 -3.0775 -2.6008 -4.4204 -5.9967 -4.6340 -4.1576 -4.6189 -4.0925 -4.3689 -5.7246 -4.3695 -3.1981 -3.9748 -4.1278 -3.7345 -3.1610 -5.0923 -4.1742 -3.9194 -4.9551 -6.6970 -4.1372 -5.7243 -4.4820 -4.8218 -3.9935

figure
plot(cgxs.', cgft.')
% plot(cgxs(:,[1 2]), cgft(:,[1 2]))
% plot(cgxs([1 2],:).', cgft([1 2],:).')
hold on
plot(cgxs50, val50, 'sr')
hold off
grid
xlabel('cgxs')
ylabel('cgft')

cgxs50 = 1×53

-4.6388 -5.5379 -5.4702 -3.4883 -3.0775 -2.6008 -4.4204 -5.9967 -4.6340 -4.1576 -4.6189 -4.0925 -4.3689 -5.7246 -4.3695 -3.1981 -3.9748 -4.1278 -3.7345 -3.1610 -5.0923 -4.1742 -3.9194 -4.9551 -6.6970 -4.1372 -5.7243 -4.4820 -4.8218 -3.9935

If you want to strictly interpolate on the ‘global’ y-value of 0.5, replace that with ‘val50’ in my code. My code scales the 50% value with the ranges of the individual ‘cgft’ vectors.

.

[Phillip]

Thank you very much!!!

(I would have never thought of that :))

[Star Strider]

As always, my pleasure!

I’m sure you would have, just as I did, when the need arose. I developed that approach when I neded to interpolate a quasi-periodic function and had to define individual interpolatioon regions for each time the signal crossed a specific threshold. It very nicely generalises to other sorts of problems, and actually seems to make the code more efficient because it narrows the interpolation region.

Answer 2

[Walter Roberson]

cgft(nsub, :) is all the same value so the second output of unique is only a single entry

[Phillip]

HI Walter - thank you for your reply. I've uploaded an updated cgft file. It will run up to subject 8 (but the output for yc is wrong from subject 3). You wouldn't know what else might be wrong?