I'm currently working on a wide gamut light source using red, green and blue LED emitters. From an internal xyY (or CIE XYZ) representation, I can reach any color or color temperature via a 3x3 transformation matrix. The matrix is calculated from the chromaticity coordinates and the relative luminance of the actual RGB emitters. This works well both in theory and in practice. 

However, the RGB LEDs emitts a discontinuous spectrum with very little energy between red and green. I want to add an orange or amber LED to improve the spectrum and the color rendering index.
Starting from CIE xyY, how do I calculate RGB plus Amber? The entire CIE model is based on tristimulus and I can't see how I can use it to calculate a fourth color.

The ideal would probably be a model that could accept any number of colors.

---

ADDED based on discussion:

![CIE xy diagram with RGB emitters](https://i.stack.imgur.com/ZfEQd.png)

Above is an illustration of how I imagine the RGB math works. 

I measure the chromaticity coordinates (x,y) of each of the Red, Green and Blue emitters, and their relative brightness (Y).

From xyY I calculate CIE XYZ. This is needed because xy plus Y is a projection, XYZ is the actual 3-dimensional color space. I put the XYZ values for the three emitters into a matrix and calculate the inverse of that matrix. This inverse matrix represents the difference between the "actual" (human vision) and the properties of these particular emitters.
If I want to display a particular color, say the white shown above to the right, I can take the desired coordinates, send them through the matrix, and get the required strength for each emitter (shown as arrows).

Now RGBA:
![CIE xy diagram with RGBA emitters](https://i.stack.imgur.com/6cQ7z.png)

Originally I thought that the math for one more emitter (RGBA) was similar to the math for RGB. That I could use a 3\*4 matrix to calculate RGBA, the same way I use a 3\*3 matrix to calculate RGB. However, it seems like there are problems calculating the inverse of a non-square matrix. Some sources say it can not be done, some say if can be done, but the resulting matrix be lacking in some respect. This is WAY beyond my understanding! :-)

@percusse suggests that a 3\*4 matrix can be used. If so, how can I calculate the inverse matrix (I'm on a shoestring budget, software like MATLAB is out of reach).

---

Second EDIT:

Based on the input from @joriki and @percusse I've tried to solve this on paper. I've spent *a lot* of paper, but I can't seem to find a way to do this that can be implemented as computer code, - or even produces the right answer! I'm probably making mistakes in the practical solving, but that is not actually critical. Computing will have to be done by a C implementation of a solving algorithm (gaussian elimination?) that is known to be good.

This would be typ XYZ values for the four emitters calculated form measured xyY coordinates (actual values will have better precision):

$$
\left[ \begin{array}{ccc}
0.47 & 0.11 & 0.19 & 0.34\\
0.20 & 0.43 & 0.11 & 0.26\\
0.00 & 0.06 & 1.12 & 0.01\end{array} \right]

\left[ \begin{array}{ccc}
R \\\
G \\\
B \\\
A \end{array} \right]
=================================================================

\left[ \begin{array}{ccc}
X \\\
Y \\\
Z \end{array} \right]
$$

I've been thinking about optimization and there are a number of parameters that affects the optimal mix, mainly spectrum, efficacy, and heat. For a small system, it is probably enough to worry about the extremes of the amber emitter (avoid max amber when emitting orange-ish light, avoid min amber when emitting any other color). A solution is already suggested by @joriki ["This selects the solution that covers the spectrum most evenly"] but I don't understand the math :-)

So I need to get this system of equations into a form that generates a single answer within 500us of computing time on a small embedded processor :-) Any guidance on how to get a step closer a practical implementation would be greatly appreciated!

---

Third EDIT:
I've set up a test that can drive 4 emitters, and a spectrometer to measure the output. The relative intensities of the emitters are tweaked to give a correlated color temperature of roughly 6000 Kelvin (midday daylight). 

RGB at ca. 6000K:

![RGB at ca. 6000K](https://i.stack.imgur.com/DIsN9.png)

RGB + Amber at ca. 6000K:

![RGB + Amber at ca. 6000K](https://i.stack.imgur.com/q8cc7.png)

RGB + White at ca. 6000K:

![RGB + White at ca. 6000K](https://i.stack.imgur.com/4JV7Z.png)

The first image shows the spectrum from 3 emitters, Red, Green and Blue. There is very little light between 560 and 610 nm. 
The next image shows the spectrum when Amber is added to RGB. Amber improves the situation significantly. (Yellow might be better, but suitable high brightness yellow LEDs can't be found).
The last image shows the spectrum when White is added to RGB. White LEDs are actually phosphor converted blue. The phosphor can be made to retransmitt over a fairly broad spectrum. This seems to give the best result in terms of even spectrum.

I *think* I have working code for Gaussian Elimination. The question now is how do I add mean square minimization to the equations in such a way that I end up with a single answer? I probably need some hints on how to solve this in practice. Sorry! :-) 

---

Fourth and fifth EDIT:

So I have measured the spectra from 380 to 780nm with 1nm resolution. The output is measured at equal input value.
![Equal input R, G, B, A and W emitters](https://i.stack.imgur.com/UAxh7.png)

I calculated the area under the curve by trigonometry. I calculated the average size for the 400 trapeziods between 380 and 780nm for R, G, B and A (values are scaled to me more manageable):

$\langle R\rangle = 19.8719507$  

$\langle G\rangle = 13.39000051$  

$\langle B\rangle = 29.30636046$  

$\langle A\rangle = 8.165754589$  

And also the average for the product of all six (plus four) combinations of emitter pairs.
I then took a stab at assembling this into a covariance matrix:

$$
\left[ \begin{array}{cccc}
43.74282392 & -2.642812728 & -5.823745503 & -0.26554119\\\
-2.642812728 & 8.563382072 & -0.969894212 & -0.946563019\\\
-5.823745503 & -0.969894212 & 62.81754221 & -2.393057209\\\
-0.26554119 & -0.946563019 & -2.393057209 & 8.136438369\end{array} \right]
$$

The matrix is assembled like this:
$$
\left[ \begin{array}{ccc}
\langle RR \rangle - \langle R \rangle\langle R \rangle & \langle RG \rangle - \langle R \rangle\langle G \rangle & \langle RB \rangle - \langle R \rangle\langle B \rangle & \langle RA \rangle - \langle R \rangle\langle A \rangle\\\
\langle RG \rangle - \langle R \rangle\langle G \rangle & ... & ... & ...\\\
... & ... & ... & ...\\\
... & ... & ... & ...\end{array} \right]
$$

Here are the measured color coordinates of the RGBA emitters, and sample values for XYZ:

$$
\left[ \begin{array}{cccc|c}
0.490449254 & 0.100440581 & 0.221653947 & 0.343906601 & 0.75\\\
0.204678363 & 0.421052632 & 0.16374269 & 0.210526316 & 1.00\\\
-0.011955512 & 0.07388664 & 1.464803251 & -0.012677086 & 0.75\end{array} \right]
$$

I've tried to get the above matrix into echelon form(?) by gaussian elimination, and then get the RGB values on the form $u + Av$ by substitution.

$R: 0.97921341 + A \* -0.701207308$   

$G: 1.730718699 + A \* -0.1767291$   

$B: 0.43 + A \* 0.012215723$   

The next step seems to be to calculate $Q$. This has been answered by @joriki, but am not used to the notation and I'm not at all sure how to translate the greek shorthand to a form where I can calculate the values. If this gets too basic for this forum, let me know and I'll take it offline. 

I have trouble understanding this calculation:

$$
\begin{eqnarray}
\mu
&=&
-\frac{\sum\_{\alpha,\beta}M\_{\alpha\beta}x\_\alpha y\_\beta}{\sum\_{\alpha,\beta}M\_{\alpha\beta}y\_\alpha y\_\beta}
\;.
\end{eqnarray}
$$

Not entirely sure what the $x$ and $y$ values are? A pointer to an example of what this $M\_{\alpha\beta}x\_\alpha y\_\beta$ look like in non-algebraic form would be very helpful.

---

Sixth EDIT:

So let me try to explain how I understand what needs to be done:
With a set of measured RGBA emitter color coordinates and an ZYX value (the color we want the emitters to generate) as input we calculate two values for each emitter. The values are 

$R = u\_{RED} + Av\_{RED}$   

$G = u\_{GREEN} + Av\_{GREEN}$   

$B = u\_{BLUE} + Av\_{BLUE}$   

$A = A$

The calculation involves gaussian elimination and substitution, and I have written code that performs those calculations.

The value of A should preferably be the one that, together with RGB, produces the most even spectrum. This appears to be the calculation of A:

$$
\begin{eqnarray}
A
&=&
-\frac{\sum\_{\alpha,\beta}M\_{\alpha\beta}x\_\alpha y\_\beta}{\sum\_{\alpha,\beta}M\_{\alpha\beta}y\_\alpha y\_\beta}
\;.
\end{eqnarray}
$$

One element of this equation is $M\_{\alpha\beta}$ which is a 4 \* 4 covariance matrix that we have precalculated from the emitter spectra. 

This is as far as I am right now. I don't understand from the above notation how the math works. Do I run every possible combination of emitter colors through the matrix add them all up? I have to admit I am completely lost! :-)