Use interp1 to interpolate a matrix row-wise

I am currently trying to expand some code to work with matrices and vectors instead of vectors and scalars. So the same calculations are to be done row-wise for n number of rows. How do I get interp1 to do this?

before I used something like this:

new_c = interp1(error,c,0,'linear',extrap')

It is used to find the value of c when an error approaches zero. Now I tried to just enter the matrices where each row is the same as the vector I used before and I get the error message "Index exceeds matrix dimensions".

I tried changing the zero to a vector of zeros but that did not help. I know I could solve it with a for-loop where I evaluate each row individually but I would prefer not to since I assume the matrix operation would save a lot of time.

Answer 1

[the cyclist]

Here is an example adapted from the online documentation ("doc interp1"):

x = 0:10; 
y1 = sin(x); 
y2 = 2*sin(x);
y = [y1;y2]';
xi = 0:.25:10; 
yi = interp1(x,y,xi); 

figure
plot(x,y,'o',xi,yi)

[Gustaf Lindberg]

Thank you for your reply but sadly it did not help. In my case, both the x and the y are matrices and I want to run interp1 row for row. In that example, x is a vector. I tried to transpose my y-matrix as it is done in the example but that made no difference. I still get an error saying that the index exceeds the matrix dimension.

[Jan]

Then please post a small running example, which reproduces your problem. I cannot follow your descriptions exactly enough to understand the problem, when you explain it in text form.

[the cyclist]

The documentation for interp1() is explicit in that x [the first input to interp1()] has to be a vector. I assumed that you had the same x values for each row of your y matrix, and that is what my example does.

If you do not have that, I'm not sure you can do this other than via a for loop. (A quick web search on the keywords suggests that it is not possible.)

The answer from Jan in this thread has a faster interpolation function than interp1(), if that helps: http://www.mathworks.com/matlabcentral/answers/44346

[Gustaf Lindberg]

I was hoping to avoid yet another loop because it's already quite a few nested loops and the interpolation is one of the most time consuming in the code. Guess I'll have to accept the extra loop.

I do not dare to try another method since I have some stability issues. It's actually not really an interpolation but more of an extrapolation. It's part of a CFD problem where I iterate through lots of cells lots of times so stability is important, even more important than speed. My supervisor has tried many different kinds of methods and he has found interp1 with the flags 'linear' and 'extrap' to be the most stable.

Thanks for all your help.

Answer 2

[Jan]

INTERP1 is slow and calling it repeatedly in a loop has a large overhead. But a linear interpolation can be implemented cheaper:

function Yi = myLinearInterp(X, Y, Xi)  
% X and Xi are column vectros, Y a matrix with data along the columns
[dummy, Bin] = histc(Xi, X);  %#ok<ASGLU>
H            = diff(X);       % Original step size

% Extra treatment if last element is on the boundary:
sizeY = size(Y);
if Bin(length(Bin)) >= sizeY(1)
   Bin(length(Bin)) = sizeY(1) - 1;
end

Xj = Bin + (Xi - X(Bin)) ./ H(Bin);
% Yi = ScaleTime(Y, Xj);  % FASTER MEX CALL HERE
% return;

% Interpolation parameters:
Sj = Xj - floor(Xj);
Xj = floor(Xj);

% Shift frames on boundary:
edge     = (Xj == sizeY(1));
Xj(edge) = Xj(edge) - 1;
Sj(edge) = 1;           % Was: Sj(d) + 1;

% Now interpolate:
if sizeY(2) > 1
   Sj = Sj(:, ones(1, sizeY(2)));  % Expand Sj
end
Yi = Y(Xj, :) .* (1 - Sj) + Y(Xj + 1, :) .* Sj;

The M-version is faster than INTERP1 already, but for the faster MEX interpolation: FEX: ScaleTime. Then the above code is 10 times faster than INTERP1.

Answer 3

[Thorsten]

 for i = 1:size(E, 1)

   new_c(i) = interp1(E(i, :), C(i, :), 0, 'linear', 'extrap');

 end

[Gustaf Lindberg]

That's how I solved it but it's not the answer to my question. I would like a way to do exactly that but without the for-loop.

[José-Luis]

What's wrong with the for loop? It is probably faster, and clearer, than the alternatives.

Answer 4

[José-Luis]

Without a loop, but slower:

nRows = size(E,1);
your_array = cell2mat(arrayfun(@(x) {interp1(E(x, :), C(x, :), 0, 'linear', 'extrap')},(1:nRows)','uniformoutput',false);

Answer 5

[Sean de Wolski]

y = toeplitz(1:10);
interp1((1:10).',y,(1:0.5:10))

[José-Luis]

The values of x change for each row of y.

[Sean de Wolski]

Ahh.

Then just use a for-loop!

Answer 6

[Matt J]

I assume 'error' is always non-negative? If so, you're really just trying to linearly extrapolate the first 2 data points in each row, which can be done entirely without for-loops and also without INTERP1,

    e1=error(:,1);
    c1=c(:,1);
    e2=error(:,2);
    c2=c(:,2);

    slopes=(c2-c1)./(e2-e1);

    new_c = c1-slopes.*e1;

[Matt J]

Note that new_c is just the y-intercepts of the line defined by the first 2 points.