Problem and attempted workaround for using interp1() with the matlab coder

Hi all, I had wanted to replicate a heating cycle for a simulation based upon given heating rates (rate = [3 2 1] K/s) across a 1-D body at a selection of points in the domain (i.e. pts = [0 .5 1]) and allowing different stoptimes for the ramprates. I have run into a few problems when trying to implement a basic approach in matlab to an algorithm that can be used along with the rest of my simulation (which is compiled with the coder due to the need for increased speed; additionally I did not want to declare extrinsic functions because of the overhead cost). The main problems result from the lack of 'cubic' or 'pchip' interpolation options

Some sample code:

    rate = [3 2 1];
    stoptime = [20 40 100];
    ramp = [0 0 0];
    pts = [0 0.5 1];
    run1 = cell(100);
    run2 = run1;

%Matlab implementation that has the desired behavior

for i=1:max(stoptime)
  %Identify the Temperature @ pts for the current time   
  ramp = i.*((stoptime>i).*rate) + (stoptime<=i).*rate.*stoptime;
  run1(i)={interp1(pts,ramp,linspace(0,1),'linear')};%Interpolate
end

%Workaround for linear interpolation

stimes = unique([0 stoptime]);
rates = cell(size(stimes));
rateint = rates;
rates(1) = {[0 0 0]};
rateint(1) = {zeros(1,100)};

for i=2:numel(stimes)
  rates(i) = {(stoptime>=stimes(i)).*rate};
  rateint(i) = {interp1([0 0.5 1],rates{i},linspace(0,1),'linear')};%Interpolate the ramp rate for each stoptime interval
end

for i=1:max(stoptime)
  %Add contributions from the portion of the current 'stimes' interval
  %and the previous intervals
  run2(i)={(i-stimes(find(i>stimes,1,'last')))...
      *rateint{find(i<=stimes,1,'first')}...
      +sum(bsxfun(@times,((stimes<i).*diff([0 stimes]))',cell2mat(rateint')))};
end

%Plot Comparison (works for 'linear', but wrong if interp1 options set to 'cubic' or 'pchip'

for i=1:max(stoptime)    
  plot(run1{i},'b+'),hold on,plot(run2{i},'ro'),ylim([0 100]),title(num2str(i))
  pause(0.1)   
  hold off    
end

This can be run to show the profile that I desire as run1 and the attempt at replication without using interp1 at each time (due to the lack of 'cubic' or 'pchip' options).

I would appreciate any guidance on correcting the rampintervals to correctly produce the desired output with the 'cubic' option. Or if there are any FEX functions/code that compatible with the coder?

Thanks, Joe

Answer 1

[Joseph]

If anyone is intrested in a version of 'pchip' and/or 'spline' that can easily be modified for use in a MEX function (because interp1() can't use those options) than a good alternative is from the Cleve Moler's 'Numerical Computing with Matlab' ( http://www.mathworks.com/moler/index_ncm.html ).

The pchiptx and splinetx can be compiled with a slight modification to the subfunctions 'pchipslopes' and 'splineslopes'. The only modification is in the size matrix size initialization to avoid an 'out of bounds' error from the compiler.

Modification:

   n = length(h)+1;
   d = zeros(1,length(h)+1);

Thank you Cleve, even if you don't know it.