Random square wave of varying individual pulse width but fixed total duration

The problem I have at hand is as follows:

A square wave (total duration: a sec) which can be (randomly) 1 for either 33 ms or 66 ms but the sum of the 1-duration is fixed (a x b x 33ms). It's 0 elsewhere. b is a fixed parameter.

I'm referring to this existing solution but the sum part in my problem is throwing me off:

https://nl.mathworks.com/matlabcentral/answers/80755-how-to-generate-a-square-wave-with-random-normal-distributed-period-and-pulse-width

[Sam Chak]

@Deepayan Bhadra

Can you sketch your desired square wave of varying width but with fixed sum of width?

[Deepayan Bhadra]

@Sam Chak: Thanks for the clarification. I framed the title incorrectly. Corrected. Hopefully, it's clear now.

[Ishaan Mehta]

Hey Deepyan

I understand that you are trying to generate square waves with following characteristics:

• total time duration: a seconds
• for each second, the duration of 1's can be b*0.033 seconds only
• the duration of each pulse must be either 0.033s of 0.066s only

I assume there is no restriction that a given pulse must start and end in the same second.

If my understanding and assumptions are valid, here is a code to generate such random waves:

a = 4;
b = 30;

t = 0:0.001:a;
y = zeros(1, a*1000 +1);

for i = 1:a
    pulsesLeftForCurrentSecond = b;

    j = (i-1)*1000 +1;
    while j <= i*1000
        if pulsesLeftForCurrentSecond <= 0
            break;
        end
        
        % pulseSpacing should decrease with an increase in value of b
        pulseSpacing = min([0.99 (0.92 + (16 - b) * 0.01)]);

        if(rand(1) > pulseSpacing) 
            random1or2 = randi([1 2]);

            if pulsesLeftForCurrentSecond <= 1
                random1or2 = 1;
            end
          
            pulseWidth = random1or2 * 33;
            y(j:j+pulseWidth-1) = 1;
            j = j + pulseWidth -1;
            pulsesLeftForCurrentSecond = pulsesLeftForCurrentSecond - random1or2;
        end

        j = j+1;
    end

    pulsesLeftForCurrentSecond
end

y(a * 1000+2: end) = [];
plot(t,y);
ylim([-1 2]);

Please note that due to the random nature of the waves, sometimes the last pulse might end up exceeding the 'a' second duration, in such cases rest of the portion is truncated and the sum of 1-duration would be less than a*b*0.033s. But such cases are extremely rare and occur only if b > 27.

(I understand that the maximum value of b can be 31 to satisfy a*b*33ms < a*1000ms)

I hope this helps while we wait for a more reliable solution.

Ishaan