T; dg
To strictly compare your implementation with Dolby specifications, you are missing a few things,
Decoder mixing levels (not relationships) are incorrect,
the LFE channel is not processed,
The direct channels of the left and right surround channels are not phase shifted , delayed or passed through HPF ,
And the center channel does not go through Bandpass .
Right combination
The clipping problem is the result of your maxtrix mixing levels (relationships are great), for example, the central channel is on amp'd with 41.4% with 0.707 + 0.707 = 1.414 ; therefore, to maintain the right relationship, simply halve your mixing levels.
With that in mind, your maxtrix should now look like this:
Left Right Center 0.354 0.354 Left 0.5 0 Right 0 0.5 SLeft 0.436 0.245 SRight 0.245 0.436
Adding LFE
If you do not intentionally leave the LFE channel, you can decode it in the same way as the center channel, but you must apply LPF at 120 Hz (Dolby standard).
Creating surround channels is actually "surround"
To process the left and right surround channels, you need to go through the left and right channels through HPF at a frequency of 100 Hz, then apply a phase shift of 90 o (then invert the left side); mix and then add a delay (I believe Dolby does not indicate the exact amount of time, but after some testing I found that "sweetspot" is about 10 milliseconds).
(If you are going to play with a delay, I would recommend doing it between 5 ms to 12.5 ms (feel free to experiment). If the delay is too small, you get a “compressed / compressed” sound mix too long, and you end up with a terrible higher frequency echoes in the background, ideally the end result should sound as “air / open” as possible, but without any hint of echo / reverb.)
leftSamples -> HPF -> phase shift -> invert -> mix \ -> delay rightSamples -> HPF -> phase shift -> mix /
Central channel
Dolby also indicates that the center channel should pass through a bandwidth of 70 Hz to 20 kHz (after mixing).
Finally...
So, putting everything together, you should get the following:
public void DolbyProLogicII(float[] leftSamples, float[] rightSamples, int sampleRate) { var ii = Math.Min(leftSamples.Length, rightSamples.Length); var c = new float[ii]; var l = new float[ii]; var r = new float[ii]; var sl = new float[ii]; var sr = new float[ii]; var lfe = new float[ii]; // Process center channel for (var i = 0; i < ii; i++) { // Best to be as precise as possible. c[i] = (leftSamples[i] * 0.35355339059327376220042218105242f) + (rightSamples[i] * 0.35355339059327376220042218105242f); } c = LinkwitzRileyHighPass(c, sampleRate, 70); c = LinkwitzRileyLowPass(c, sampleRate, 20000); // Process front left channel for (var i = 0; i < ii; i++) { l[i] = leftSamples[i] * 0.5f; } // Process front right channel for (var i = 0; i < ii; i++) { r[i] = rightSamples[i] * 0.5f; } // Process left samples for SL channel var slL = new float[ii]; for (var i = 0; i < ii; i++) { slL[ii] = leftSamples[i] * 0.43588989435406735522369819838596f; } slL = LinkwitzRileyHighPass(slL, sampleRate, 100); slL = PhaseShift(slL, sampleRate, true); // Process right samples for SL channel. var slR = new float[ii]; for (var i = 0; i < ii; i++) { slR[i] = rightSamples[i] * 0.24494897427831780981972840747059f; } slR = LinkwitzRileyHighPass(slR, sampleRate, 100); slR = PhaseShift(slR, sampleRate); // Combine new left & right samples for SL channel for (var i = 0; i < ii - 1; i++) { sl[i] = slL[i] + slR[i]; } sl = Delay(sl, sampleRate, 10); // Process left samples for SR channel var srL = new float[ii]; for (var i = 0; i < ii; i++) { srL[i] = leftSamples[i] * 0.24494897427831780981972840747059f; } srL = LinkwitzRileyHighPass(srL, sampleRate, 100); srL = PhaseShift(srL, sampleRate, true); // Process right samples for SR channel var srR = new float[ii]; for (var i = 0; i < ii; i++) { srR[i] = rightSamples[i] * 0.43588989435406735522369819838596f; } srR = LinkwitzRileyHighPass(srR, sampleRate, 100); srR = PhaseShift(srR, sampleRate); // Combine new left & right samples for SR channel for (var i = 0; i < ii - 1; i++) { sr[i] = srL[i] + srR[i]; } sr = Delay(sr, sampleRate, 10); // Process LFE channel. for (var i = 0; i < ii; i++) { lfe[i] = (leftSamples[i] * 0.35355339059327376220042218105242f) + (rightSamples[i] * 0.35355339059327376220042218105242f); } lfe = LinkwitzRileyLowPass(lfe, sampleRate, 120); } public float[] PhaseShift(float[] samples, double sampleRate, bool invertOutput = false) { var depth = 4.0; var delay = 100.0; var rate = 0.1; var newSamples = new float[samples.Length]; double wp, min_wp, max_wp, range, coef, sweepfac; double inval, x1, outval = 0.0; double lx1 = 0.0, ly1 = 0.0, lx2 = 0.0, ly2 = 0.0, lx3 = 0.0, ly3 = 0.0, lx4 = 0.0, ly4 = 0.0; // calc params for sweeping filters wp = min_wp = (Math.PI * delay) / sampleRate; range = Math.Pow(2.0, depth); max_wp = (Math.PI * delay * range) / sampleRate; rate = Math.Pow(range, rate / (sampleRate / 2)); sweepfac = rate; for (var i = 0; i < samples.Length; i++) { coef = (1.0 - wp) / (1.0 + wp); // calc coef for current freq x1 = (inval = (double)samples[i]); ly1 = coef * (ly1 + x1) - lx1; // do 1st filter lx1 = x1; ly2 = coef * (ly2 + ly1) - lx2; // do 2nd filter lx2 = ly1; ly3 = coef * (ly3 + ly2) - lx3; // do 3rd filter lx3 = ly2; ly4 = coef * (ly4 + ly3) - lx4; // do 4th filter lx4 = ly3; // final output outval = ly4; if (invertOutput) { newSamples[i] = -(float)outval; } else { newSamples[i] = (float)outval; } wp *= sweepfac; // adjust freq of filters if (wp > max_wp) // max? { sweepfac = 1.0 / rate; // sweep back down } else { if (wp < min_wp) // min? { sweepfac = rate; // sweep back up } } } return newSamples; } public float[] Delay(float[] samples, int sampleRate, double milliseconds) { var output = new List<float>(samples); var ii = (sampleRate / 1000) * milliseconds; for (var i = 0; i < ii; i++) { output.Insert(0, 0); } return output.ToArray(); } public float[] LinkwitzRileyHighPass(float[] samples, int sampleRate, double cutoff) { if (cutoff <= 0 && cutoff >= sampleRate / 2) { throw new ArgumentOutOfRangeException("cutoff", "The cutoff frequency must be between 0 and \"sampleRate\" / 2."); } if (sampleRate <= 0) { throw new ArgumentOutOfRangeException("sampleRate", "The sample rate must be more than 0."); } if (samples == null || samples.Length == 0) { throw new ArgumentNullException("samples", "\"samples\" can not be null or empty."); } var newSamples = new float[samples.Length]; var wc = 2 * Math.PI * cutoff; var wc2 = wc * wc; var wc3 = wc2 * wc; var wc4 = wc2 * wc2; var k = wc / Math.Tan(Math.PI * cutoff / sampleRate); var k2 = k * k; var k3 = k2 * k; var k4 = k2 * k2; var sqrt2 = Math.Sqrt(2); var sq_tmp1 = sqrt2 * wc3 * k; var sq_tmp2 = sqrt2 * wc * k3; var a_tmp = 4 * wc2 * k2 + 2 * sq_tmp1 + k4 + 2 * sq_tmp2 + wc4; var b1 = (4 * (wc4 + sq_tmp1 - k4 - sq_tmp2)) / a_tmp; var b2 = (6 * wc4 - 8 * wc2 * k2 + 6 * k4) / a_tmp; var b3 = (4 * (wc4 - sq_tmp1 + sq_tmp2 - k4)) / a_tmp; var b4 = (k4 - 2 * sq_tmp1 + wc4 - 2 * sq_tmp2 + 4 * wc2 * k2) / a_tmp; var a0 = k4 / a_tmp; var a1 = -4 * k4 / a_tmp; var a2 = 6 * k4 / a_tmp; var a3 = a1; var a4 = a0; double ym1 = 0.0, ym2 = 0.0, ym3 = 0.0, ym4 = 0.0, xm1 = 0.0, xm2 = 0.0, xm3 = 0.0, xm4 = 0.0, tempy = 0.0; for (var i = 0; i < samples.Length; i++) { var tempx = samples[i]; tempy = a0 * tempx + a1 * xm1 + a2 * xm2 + a3 * xm3 + a4 * xm4 - b1 * ym1 - b2 * ym2 - b3 * ym3 - b4 * ym4; xm4 = xm3; xm3 = xm2; xm2 = xm1; xm1 = tempx; ym4 = ym3; ym3 = ym2; ym2 = ym1; ym1 = tempy; newSamples[i] = (float)tempy; } return newSamples; } public float[] LinkwitzRileyLowPass(float[] samples, int sampleRate, double cutoff) { if (cutoff <= 0 && cutoff >= sampleRate / 2) { throw new ArgumentOutOfRangeException("cutoff", "The cutoff frequency must be between 0 and \"sampleRate\" / 2."); } if (sampleRate <= 0) { throw new ArgumentOutOfRangeException("sampleRate", "The sample rate must be more than 0."); } if (samples == null || samples.Length == 0) { throw new ArgumentNullException("samples", "\"samples\" can not be null or empty."); } var newSamples = new float[samples.Length]; var wc = 2 * Math.PI * cutoff; var wc2 = wc * wc; var wc3 = wc2 * wc; var wc4 = wc2 * wc2; var k = wc / Math.Tan(Math.PI * cutoff / sampleRate); var k2 = k * k; var k3 = k2 * k; var k4 = k2 * k2; var sqrt2 = Math.Sqrt(2); var sq_tmp1 = sqrt2 * wc3 * k; var sq_tmp2 = sqrt2 * wc * k3; var a_tmp = 4 * wc2 * k2 + 2 * sq_tmp1 + k4 + 2 * sq_tmp2 + wc4; var b1 = (4 * (wc4 + sq_tmp1 - k4 - sq_tmp2)) / a_tmp; var b2 = (6 * wc4 - 8 * wc2 * k2 + 6 * k4) / a_tmp; var b3 = (4 * (wc4 - sq_tmp1 + sq_tmp2 - k4)) / a_tmp; var b4 = (k4 - 2 * sq_tmp1 + wc4 - 2 * sq_tmp2 + 4 * wc2 * k2) / a_tmp; var a0 = wc4 / a_tmp; var a1 = 4 * wc4 / a_tmp; var a2 = 6 * wc4 / a_tmp; var a3 = a1; var a4 = a0; double ym1 = 0.0, ym2 = 0.0, ym3 = 0.0, ym4 = 0.0, xm1 = 0.0, xm2 = 0.0, xm3 = 0.0, xm4 = 0.0, tempy = 0.0; for (var i = 0; i < samples.Length; i++) { var tempx = samples[i]; tempy = a0 * tempx + a1 * xm1 + a2 * xm2 + a3 * xm3 + a4 * xm4 - b1 * ym1 - b2 * ym2 - b3 * ym3 - b4 * ym4; xm4 = xm3; xm3 = xm2; xm2 = xm1; xm1 = tempx; ym4 = ym3; ym3 = ym2; ym2 = ym1; ym1 = tempy; newSamples[i] = (float)tempy; } return newSamples; }