The Mueller optical calculus is well suited to handling this kind of problem, since the LEDs are intensity modulated. Using my optical calculus software, here is a simple computer simulation that assumes three independently modulated red, green and blue LEDs, three separate bandpass optical filters and one photodiode followed by a transimpedance amplifier:
This is my screenshot of my simulation program. It is not a drawing.
In this simulation, I assume that the wavelength space is discrete, consisting of 401 colors: 400 to 800 nm, in increments of 1 nm. The simulation runs from t = 0 to t = 1.6383 s, in increments of 0.0001 s. (This was to facilitate doing a 16k DFFT at the end of the simulation.) I used the following (not so realistic) LED spectra in the LED "light bulbs" in the simulation:
This is my screenshot of the LED spectra. The LED intensities are arbitrarily scaled and are purely synthetic: I did not have real digitized spectra available. These could trivially be replaced with real LED spectra: it is just an import of ASCII text files. The photodiode responsivity is shown by the black trace, with scale on the right. The units are A/W.
For the bandpass optical filters, I synthesized the following deliberately overlapping transmittance spectra:
This is my screenshot of the filter bandpass spectra and the photodiode's responsivity. For convenience, the photodiode responsivity is repeated: shown by the black trace, with scale on the right.
Only a small amount of noise was added: 10 \$\mu \$V electrical noise after the preamp and 0.01% optical noise on each of the three color channels. The transimpedance was 10 k\$\Omega \$. There is no difficulty in adding noises. Rather, the difficulties are in knowing what noises to add, how much of each one to add, and where to add them. This is not trivial.
The two mirrors and two beamsplitting cubes were specified as being ideal, though that is not required. The three LED modulators are simply producing unipolar squarewaves at 100, 175 and 375 Hz.
When the simulation runs, which takes just over one minute, the voltage output is as shown below:
This is my screenshot of the temporal output, from the little oscilloscope. You get this when you double-click on the oscilloscope block.
The PSD is computed at the end of the single simulation and is shown below on a linear-linear scale:
This is my screenshot of the power spectral density output, from the little PSD block. You get this when you double-click on the PSD block.
The three fundamental frequencies, i.e., 100, 175 and 375 Hz, are clearly present, along with several squarewave harmonics.
In the above simulation, it was assumed that the LED emissions were unpolarized, but the full Mueller calculus was used anyway because there was nothing to be gained by simplification of the Mueller matrices and Stokes vectors.
So where is the math? Actually, it is all done in the blocks. An optical calculus is really just a little bit of linear algebra. For any one color, a Stokes 4-vector arrives at the input to a block, gets left-multiplied by the block's resident 4x4 Mueller matrix, and then the resulting Stokes 4-vector is sent out. For any given time step in the simulation, all Stokes vectors arriving at a given block are processed at once: an arriving stack of Stokes vectors is left-multiplied by a resident stack of Mueller matrices, producing an outgoing stack of Stokes vectors. See the Addendum below for the stack exchange format for the optical calculus blocks.
All this is intentionally invisible to the user, unless they happen to want to look at the commented source code. In that case, they would place the cursor on the block, hold down the option key, and double-click. Then the full commented source code, and all the rest, would be available for inspection, modification, etc.
However, if the light is unpolarized, then simplification is very simple: the situation is exactly as alex.forencich stated in his answer. The special case equations are then like this:
These are the equations I constructed using MathType and annotated with Keynote (all on a Mac). Here, i indexes the N+1 discrete colors, j indexes the evenly spaced simulation times and k indexes the three color channels. I do not know Python, but maybe this all will help a bit when you get to the programming.
The software I used consists of a commercial simulation program (Extend 6.0.8, from Imagine That, Inc.) together with my own free libraries of add-on blocks (named LightStone, my punning play on "optical calculus"). Unfortunately, as I found out fairly recently, my free libraries of blocks, developed and evolved since 1990, no longer work with the current ExtendSim program (versions 10 through 10.0.6). But they work fine with older versions of ExtendSim, e.g., ExtendSim 9.2 on a Window 10 PC. I hope to get this resolved because I have 30 years riding on it.
Update: My LightStone libraries now work with ExtendSim CP version 10.0.8, which is the latest version. They also work with all older versions of ExtendSim. Glad to have this resolved!
Addendum: The format for optical block stack exchanges, i.e., outputs and inputs of stacks of optical vectors, is shown in the next figure, which I drew about 25 years ago using (long dead now) Canvas 5.0.3:
My avatar happens to be the right hand side of Fig. 7. In this answer, there are 401 rows, i.e., discrete wavelengths. Since Mueller calculus is used, only the first column is used, i.e., the Stokes vector. Since the light is unpolarized, only the top element in the Stokes vector, i.e., "a" in the figure, is used. Of course, this is the light intensity element and it is both time and wavelength dependent: in Fig. 6, it is labeled as "Modulated LED intensity".
So it is stacks of optical vectors that are passed from optical block to optical block, kind of like this little animated cartoon GIF I made many years ago:
In actuality, the cartoon is a little misleading because stacks pass completely with each simulation time step. But it makes the important point that three dimensional arrays are passing, not just single real numbers.
