Reconstructing a Random Step Sequence using FFT

Dear All,

It may seem like a trivial question, but can't seem to figure it out.

I am using fft to reconstruct a randomly generated step sequence (0's and 1's) through Fourier series coefficients. All else is fine except my reconstructed signal is off by one sample, meaning it is shifted back by 1 sample. Can you please tell me what could I be doing wrong here?

Thank you for your time and help.

Regards, Adeel

P.S. you can skip to part 2 in the code where I am calculating the Fourier series coefficients.

%%Part 1
%%To obtain fourier series coefficients of a peroidic signal using fft
clear all; close all; clc;
dur = 0.5;    % Signal duration (sec)
step = 1e-3;  % step size
f = 50;       % frequency of the signal (Hz)
t = 0:step:dur-step;
M = length(t);
y1 = zeros(1,M);
p = 1;
%%This for loop generates a random step sequence with different 
% duty cycle each time it runs. It emulates the behavior of 
% sinusoidal PWM for use in DC to AC converters.
for i = 1:M
    n = randi([0,1],1);     % generate a random 0 or 1
    k = randi([10,50],1);   % generate a random integer between 10 & 50
    if n == 0;          % if n is 0, then the k entries of y1 will be 0
        for j = 1:k
            y1(p) = 1;
            p = p+1;    % increment index of y1 vector
        end
    else                % if n is 1, then the k entries of y1 will be 1
        for j = 1:k
            y1(p) = 0;
            p = p+1;    % increment index of y1
        end

      end
      if p >=M
          y1 = y1(1:M);
          break;
      end
  end
  figure(1)
  plot(t,y1,'LineWidth',2)
  axis([0 dur -0.1 1.1])

%%----------------------------------------------
%%Part 2
%%Calculate the Fourier series coefficients of y1
N = length(y1);     % window length/data length
Y1 = fft(y1)/N;     % take fft
N1 = ceil(N/2);     
x = 1:1:N;
k = 0:N-1;
omega = k*2*pi*f;   % fundamental frequencies and its harmonics

an = 2*real(Y1);
an(1) = an(1)/2;
bn = -2*imag(Y1);

fapp = an(1)*ones(size(y1));
L = N1;
for k = 1:L
    fapp = fapp + an(k+1)*cos(2*pi*x*k/N) + bn(k+1)*sin(2*pi*x*k/N);
end

plot(t,y1,'-.b'); hold on
plot(t,fapp,'-r')

Answer 1

[Adeel]

Resolved. I changed x from 0:1:N-1 instead of 1:1:N, and it fixed the issue.