Pulstran function to generate biphasic pulses

How can I correct following code to get pulses like on the picture? I don't get correct shape of pulses.

Parameters: N = 400 (number of pulses), Tp = 1 (pulse width), d1 = 1 (interphase delay [us]), d2 = 1 (interpulse delay [us]). Let's say the U0=1V. It should look like this:

number_of_points_on_interval = 1000;

t = linspace(0, 1.6, 10000);

Tp = 0.001; %pulse width [ms]
d1 = 0.001; %interphase delay [ms]
d2 = 0.001; %interpulse delay [ms]

dPos = 0.001:(d1+Tp+d2):1.6; % times when positive pulses start

dNeg = dPos + Tp + 0.001; % times when negative pulses start

dT = [dPos dNeg]'; % times of pulses

% Amplitude of positive and negative pulses
U=1;
dAmp = [U*ones(numel(dPos),1); -1 * U*ones(numel(dNeg),1)];

% generating a signal with rectangular pulses
y = pulstran(t,[dT dAmp],'rectpuls',Tp);
plot(t,y)
ylim padded

Answer 1

[Mathieu NOE]

hello

have to admit I almost never used pulstran yet , so I simply prefer to come back to a very simple solution with repmat if that suffice for the job

NB that to test the code I reduced the final time for the plot - simply put back t_final = 1.6 on your side

also I changed interpulse from 1 to 3 ms just to check the code - here also change back to your own settings

t_final = 1.6/50; % <= put back t_final = 1.6 on your side 
dt = 1e-4;
t = (0:dt:t_final);

Tp = 0.001; %pulse width [ms]
d1 = 0.001; %interphase delay [ms]
d2 = 0.003; %interpulse delay [ms]

% Amplitude of positive and negative pulses
U=1;

% d2 period (zero amplitude)
t1 = (0:dt:d2-dt);
y1 = zeros(size(t1));

% pos pulse (+1 amplitude)
t2 = t1(end) + (dt:dt:Tp);
y2 = U*ones(size(t2));

% interphase (zero amplitude)
t3 = t2(end) + (dt:dt:d1);
y3 = zeros(size(t3));

% neg pulse (-1 amplitude)
t4 = t3(end) + (dt:dt:Tp);
y4 = -U*ones(size(t4));

dT = [t1 t2 t3 t4]'; % times of pulses
dAmp = [y1 y2 y3 y4]';

% how many cycles do we need to repeat the pattern ? 
repet = fix(t_final/t4(end));
y = repmat(dAmp,repet,1);
t = dt*(0:numel(y)-1);

% generating a signal with rectangular pulses
plot(t,y)

% ylim padded

[Marina Babic]

Hello! Thank you so much for your answer, but sadly this code is not helping me much.. The reason I am using pulstran is because I need a code where time variable is a part of variable "y" (in this case), because I need to incorporate time from ode15s function into the code above.

[Mathieu NOE]

well , you have t and y in my solution....see the end of my code

y = repmat(dAmp,repet,1);

t = dt*(0:numel(y)-1);

[Marina Babic]

Yes, I have seen that. But I need "t" to be a part of the y.. For instance:

bipolar_pulses = Iext*square(2*pi*frequency*t)
pulses = pulstran(t,d,x);
x = gauspuls(t,10e3,0.5); 

In all of the above, variable t is a part of the function that generates pulses. Another function will give me variable t and that t I will need to include in pulse code. I hope I was clear in explaining..

[Mathieu NOE]

ok - maybe it's because it's monday .... I was slow to understand your comment

now , you simply pass t as argument to my newly created function and there you go ...

t_final = 1.6/40; 
dt = 1e-4;
t = (0:dt:t_final);

% generating a signal with rectangular pulses
y = mypulsetrain(t);

plot(t,y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = mypulsetrain(t)

    dt = mean(diff(t));
    t_final = t(end);
    
    Tp = 0.001; %pulse width [ms]
    d1 = 0.001; %interphase delay [ms]
    d2 = 0.003; %interpulse delay [ms]

    % Amplitude of positive and negative pulses
    U=1;

    % d2 period (zero amplitude)
    t1 = (0:dt:d2-dt);
    y1 = zeros(size(t1));

    % pos pulse (+1 amplitude)
    t2 = t1(end) + (dt:dt:Tp);
    y2 = U*ones(size(t2));

    % interphase (zero amplitude)
    t3 = t2(end) + (dt:dt:d1);
    y3 = zeros(size(t3));

    % neg pulse (-1 amplitude)
    t4 = t3(end) + (dt:dt:Tp);
    y4 = -U*ones(size(t4));

    dT = [t1 t2 t3 t4]'; % times of pulses
    dAmp = [y1 y2 y3 y4]';

    % how many cycles do we need to repeat the pattern ? 
    repet = fix(t_final/t4(end));
    y = repmat(dAmp,repet,1);
    % padd y with zeros to match t length
    p = numel(t) - numel(y);
    y = [y;zeros(p,1)];


end

[Marina Babic]

Thank you for your effort! Sadly I am not able to use this code successfully with my alredy written (second part of) code:( on it's own, it works but when combined with other part of code, it gives me an error for the repmat. I've tried to fix but error isn't going. Anyway I very much appreciate the time and help from you!