top of page

JMO_Spectrum: A Matlab spectrum library

Writer's picture: Julius MuschaweckJulius Muschaweck

Updated: Mar 29, 2021

We deal routinely with spectra in illumination optics. We need to analyze LED spectra, compute color coordinates and color rendering values from spectra, integrate them, add spectra, multiply spectra with scalar weights and with other spectra (like transmission spectra), interpolate a given measured spectrum with non-equidistant wavelength values to a regular 1 nm array, and so on. In practice, this is tedious: spectra come in various formats, and the problem of dealing with two spectra als tabulated values which have two different sets of wavelengths is annoying. And while the various algorithms to compute things with spectra should not be rocket science, their implementation is sometimes tricky. Over the years, I settled on a series of Matlab routines I wrote to do all these chores for me, nice and easy. And now, I'm making the code freely available. So, this open source Matlab library is designed to make dealing with spectra easy and transparent. It is partly compatible with GNU Octave, and should be fully compatible with Matlab for Mac and Linux. Available for free download on https://github.com/JuliusMuschaweck/JMO_Spectrum.


The full documentation (see the ReadMe) contains all the details, but here is the flavor:

In this library, a spectrum is a struct which

  • contains a field named lam, which is a 1D real array with length >= 2, and strictly ascending values,

  • contains a field named val, which is a 1D real array with same length as lam.

That's it. Simple to create, by one of the two following ways:


% // a flat spectrum, known as CIE standard illuminant E
clear s1;
s1.lam = [360 830];
s1.val = [1 1];
% // the same spectrum, created with a convenience helper function that 
% // includes sanity checks
s2 = MakeSpectrum( [360 830], [1 1]); % // the same spectrum

I decided against making my spectra a Matlab class, because I wanted to make it simple to add fields in a flexible and extensible way:


s.name = "CIE standard illuminant E";
s.hopp = "topp"; % // anything you want
s.CCT = CCT(s); % // a common idiom: compute something, add it as field

This requires a tiny bit of discipline on the user's side, but this is outweighed by the simplicity, I believe.

There are routines to add and multiply spectra, where the "interweaving" of non-equal wavelength arrays is taken care of by the library.

You can create standard spectra like Planck blackbody, all the CIE standard illuminants, (bell shaped) Gauss spectra, typical white LED spectra and more. You can read and write LightTools® spectrum files, which are supplied by LED vendors, and also simple tabulated ASCII files.

And you can do a lot of colorimetric calculations, like CIE 1931 XYZ color coordinates, correlated color temperature (CCT), dominant wavelength, color rendering index, and more.

All colorimetric calculations are based on official CIE data, now available from https://cie.co.at/data-tables . Those official spectrum tables used to be buried in annexes of expensive, not freely available CIE standards; I commend CIE for making these datasets now freely available to the community, as .csv files with associated .json metadata. I imported the CIE data sets with metadata for easy access from within Matlab, ensuring data integrity using CIE-provided MD5 checksums, sample rows and sums of columns.

I did my very best to get some tricky implementation details right. The Planck locus, for example, is computed as tabular values using a fine grid of absolute temperatures, the original formula by Planck and CODATA2018 values for the physical constants. Then, I provide a spline interpolation function for the color coordinates, which is then used to find the CCT of a spectrum with a root finding algorithm with high precision:

lam = 360:830;
T = 5761.68;
ps = PlanckSpectrum(lam, T);
CCT(ps)-T

returns 0.00017 K.

Or, take linear interpolation, which is at the core of the "interweaving": For Matlab on Windows, I created a C++ DLL which is much faster than Matlab's interp1 function, to which the code falls back on other platforms.

A picture says more than a thousand words... I provide a number of plotting routines, e.g., the familiar CIE x-y horseshoe diagram:

PlotCIEBorder(ColorFill=true);

But this is much more than a graphic: It's an active Matlab figure, to which you can add your own elements, like the sRGB color gamut:

fh = figure(); clf; hold on; 
[ah, fh] = PlotCIExyBorder('Figure', fh, 'LineSpec','m',...
    'PlotOptions',{'LineWidth',1.5},'Ticks', [510 520 530],...
    'TickFontSize',10, 'ColorFill',true);
grid on;
axis equal;
axis([-0.05,0.8,-0.02,0.9]);
R = sRGB_to_XYZ(1,0,0); 
G = sRGB_to_XYZ(0,1,0); 
B = sRGB_to_XYZ(0,0,1);
gamut_x = [R.x, G.x, B.x, R.x];
gamut_y = [R.y, G.y, B.y, R.y];
plot(gamut_x, gamut_y,Color='k',Marker='o',LineWidth=2);
xlabel('CIE x'); ylabel('CIE y');
title('CIE 1931 diagram with sRGB gamut');

There is much more, like TLCI, CIEDE2000, IES TM30, CIELAB / CIELUV ...

Enjoy using the code, and please send me feedback to julius@jmoptics.de when you do.

611 views0 comments

Comentarios


bottom of page