title: Spectral data to XYZ color model in MATLAB
body: <p>I'm working in MATLAB, trying to take the reflection spectrum of a sample and display the perceived color, as an approximation.</p>
  
  <p>I used the 2-deg XYZ CMFs transformed from the CIE (2006) 2-deg LMS cone fundamentals with 0.1nm step size from <a href="http://cvrl.ioo.ucl.ac.uk/cmfs.htm" rel="nofollow noreferrer">http://cvrl.ioo.ucl.ac.uk/cmfs.htm</a> as my color matching functions.</p>
  
  <p>I used the three formulae shown here: <a href="http://www.fourmilab.ch/documents/specrend/" rel="nofollow noreferrer">http://www.fourmilab.ch/documents/specrend/</a></p>
  
  <p>I think it should be pretty straightforward, but somehow I keep getting Z values around 120 for the spectrum I'm evaluating, which is higher than the allowed range for Z.</p>
  
  <p>To get higher accuracy, I tried interpolating my matching functions, but I got very similar results than what I did previously.</p>
  
  <pre><code>function [X,Y,Z, RGB] = spectrumtoXYZ(P, Lambda, x, y, z)
  
  % - P is the spectrum
  % - Lambda is the wavelength for which I have my matching function data
  %   points
  % - x,y,z are the matching functions
  
  X=0;Y=0;Z=0;
  
  %Delta Lambda as specified here https://www.fourmilab.ch/documents/specrend/
  del = (P(end,1)-P(1,1))/length(P);
  
  for i=1:length(P)   
      %I look for the closest wavelength available in the matching function
      [~, j] = min(abs(P(i,1)-Lambda));
      fx=0;fy=0;fz=0;
  
      % Linear interpolation of Lambda for better integral accuracy
      if (P(i,1)-Lambda(j))&gt;0 &amp;&amp; j~=length(Lambda)
          fx=x(j)+(x(j+1)-x(j))/(Lambda(j+1)-Lambda(j))*(P(i,1)-Lambda(j));
          fy=y(j)+(y(j+1)-y(j))/(Lambda(j+1)-Lambda(j))*(P(i,1)-Lambda(j));
          fz=z(j)+(z(j+1)-z(j))/(Lambda(j+1)-Lambda(j))*(P(i,1)-Lambda(j));
      end
  
      if (P(i,1)-Lambda(j))&lt;=0 &amp;&amp; j~=1
          fx=x(j)+(x(j)-x(j-1))/(Lambda(j)-Lambda(j-1))*(P(i,1)-Lambda(j));
          fy=y(j)+(y(j)-y(j-1))/(Lambda(j)-Lambda(j-1))*(P(i,1)-Lambda(j));
          fz=z(j)+(z(j)-z(j-1))/(Lambda(j)-Lambda(j-1))*(P(i,1)-Lambda(j));
      end
  
      %Integral
      X=del*fx*P(i,2)+X;
      Y=del*fy*P(i,2)+Y;
      Z=del*fz*P(i,2)+Z;
  end
  
  % X,Y,Z to RGB; Ref. http://en.wikipedia.org/wiki/CIE_1931_color_space
  M=[0.41847 -0.15866 -0.082835; -0.091169 0.25243 0.015708; 0.00092090 -0.0025498 0.17860];
  RGB = M*[X;Y;Z];
  end
  </code></pre>
  
  <p>Is there a blatant mistake in my reasoning/code? I've been staring at this for hours. Any ideas of what I could test to troubleshoot?</p>