loop error with "unique" and "interp1"

1 次查看(过去 30 天)
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.)
  3 个评论
Torsten
Torsten 2023-4-30
Both arrays must be of the same size, but they aren't (see above).
Phillip
Phillip 2023-4-30
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).

请先登录,再进行评论。

采纳的回答

Star Strider
Star Strider 2023-4-30
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.
.
  2 个评论
Phillip
Phillip 2023-4-30
编辑:Phillip 2023-4-30
Thank you very much!!!
(I would have never thought of that :))
Star Strider
Star Strider 2023-4-30
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.

请先登录,再进行评论。

更多回答(1 个)

Walter Roberson
Walter Roberson 2023-4-30
cgft(nsub, :) is all the same value so the second output of unique is only a single entry
  1 个评论
Phillip
Phillip 2023-4-30
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?

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 NaNs 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by