# Numerical error while implementing a periodic pulse signal

I'm trying to illustrate the principle of constructing a periodic signal using an elementary pattern. This is the code I use (python 2.7 + numpy + scipy.signal) :

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as sfft
import scipy.signal as spsig

plt.close('all')

npts = 2**14 #number of points for both time-domain signals and FFT's
tsim = 100e-9 #ending time (in s)
fs = npts/tsim #sampling frequency

t = np.linspace(0,tsim,npts) #time bins vector
t1 = 30e-12 #pulse rise time
t2 = 80e-12 #pulse fall time
delta_t = 200e-12 #pulse width (including t1 and t2)
fsig = 1e9 #pulse repetition frequency
A = 1. #pulse amplitude

# here I construct the elementary pulse signal, it works fine
sig_elem = np.zeros(t.size)
sig_elem[(t>=0)&(t<t1)] = (A/t1)*t[(t>=0)&(t<t1)]
sig_elem[(t>=t1)&(t<delta_t-t2)] = A
sig_elem[(t>=delta_t-t2)&(t<delta_t)] = (-A/t2)*t[(t>=delta_t-t2)&(t<delta_t)] + (A/t2)*delta_t

periode = 1/fsig #this is the pulse repetition period
comb = np.reshape(np.array([np.round(t*1e9, 2)%1==0]), (npts,)) #this is where I construct the comb, I use numerical tricks to make it work
sig = spsig.fftconvolve(sig_elem, comb, 'full') #and the convolution itself, which should give a periodization of my elementary signal

##### FIGURE ####
plt.figure(1)
plt.plot(t*1e9, sig_elem, 'r',lw = 0.8)
plt.xlabel('t (ns)')
plt.xlim(-10,100)
plt.ylim(-1,5)
plt.grid(True)

# single sided FFT of the elementary signal. I restrict the plot frequency at 50 GHz since the sampling frequency is insanely high
sig_elem_fft = np.abs(sfft.fftshift(sfft.fft(sig_elem, npts)))
plt.figure(2)
plt.plot(sfft.fftshift(sfft.fftfreq(npts, 1/fs))*1e-9, sig_elem_fft, 'b',lw = 0.8)
plt.xlabel('f (GHz)')
plt.xlim(0,50)
plt.grid(True)

#plot of the dirac "comb". It seems good
plt.figure(3)
plt.scatter(t*1e9, comb)
plt.grid(True)
plt.xlim(-10,110)
plt.ylim(-0.5,1.5)

# plot of the convolve result, this is where I don't understand what's happening
plt.figure(4)
plt.plot(sig, 'r',lw = 0.8)
plt.xlabel('t (ns)')
#plt.xlim(-10,100)
#plt.ylim(-1,5)
plt.grid(True)


What I'm doing is defining a trapezoidal pulse :

$$\begin{cases} s(t) = \frac{At}{t_1} & \text{ if }0 \leq t < t_1\\ A & \text{ if } t_1\leq t < \Delta t-t_2\\ \frac{A}{t_2}\left(\Delta t - t\right) & \text{ if } \Delta t-t_2 \leq t < \Delta t \end{cases}$$

$t_1$, $t_2$ and $\Delta t$ being respectively rise time, fall time and total pulse width. A is an arbitrary amplitude. I implement it on $2^{14}$ (16384) points, making the total duration equal to 100 ns. Which makes the sampling frequency around 164 GHz.

This is what I get, the rest of the signal being equal to zero.

After creating this elementary pulse, I create a "dirac comb", using the same number of points, and using a 3 GHz repetition frequency. Then, both signals are convolved using an FFT method (meaning both signals FFT's are calculated, multiplied element-wise and the result IFFT is computed afterwards).

What I should get is a periodic signal consisting of $s(t)$ defined above being repeated each 1 ns. However, I get a strange signal, where the amplitude is sometimes A, sometimes 2*A. The period is however consistent.

Your "numeric trick" is not working. After running above code, have a look at :

plt.figure(); plt.plot(comb); plt.xlim((324,336)); plt.ylim((-0.1, 1.1))


Your "Dirac" train has peaks with width 2, hence you have the double-height amplitude pulses.

You can fix this by e.g. doing

comb2 = np.zeros(npts)
step = 1e-9*fs  # distance between peaks in terms of samples (floating point)
peakPos = step*np.arange(int(npts / step))  # integer position of the peaks
comb2[np.round(peakPos).astype(int)] = 1  # set the peaks
plt.figure()
plt.plot(comb)
plt.plot(comb2); plt.ylim((-0.1, 1.1)); plt.xlim((324, 336))