pythonsignalsfftsimulation

Incorrect modeling of the spectrum of a sequence of pulse signals

I have a task to simulate a sequence of rectangular pulses. As a result, I wrote Python code, but the end results are not encouraging.

#!/bin/python

import numpy as np
from numpy.fft import fft
import matplotlib.pyplot as plt

A = np.random.randint(1, 6) # random amplitude (from 1 to 5)
tilda = 0.5 # pulse duration
SR = 1000 # sampling rate
t = np.linspace(0, 3, SR) # time range from 0 to 2 seconds

n_pulses = 3 # number of pulses
t_i = 0.5 # interval between pulses
pulse1_start = 0.25 # the beginning of the first pulse

rect_pulses = np.zeros_like(t)
for i in range(n_pulses):
    pulse_start = pulse1_start + i * (tilda + t_i)
    pulse_end = pulse_start + tilda
    rect_pulses = np.where((t >= pulse_start) & (t <= pulse_end), A, rect_pulses)

FFT = fft(rect_pulses) # magnitude of FFT
FFT_n = FFT / len(FFT)  # normalization
f_ax = np.linspace(0, (SR), len(FFT)) # frequency axis 
half_fft = len(FFT) // 2 # half the length of the FFT
f_half = f_ax[:half_fft] # frequencies up to f_ax / 2
FFT_half = np.abs(FFT_n[:half_fft]) # FFT module up to f_ax/2
FFT_half *= 2 # doubling the magnitude 

plt.figure(figsize=(10, 10))
plt.subplot(2, 1, 1)
plt.plot(t, rect_pulses, linewidth=1.7)
plt.title('Rect. Pulse', fontsize=21)
plt.xlabel('t', fontsize=21, fontweight='bold')
plt.ylabel('x(t)', fontsize=21, fontweight='bold')
plt.grid()
plt.subplot(2, 1, 2)
plt.plot(f_half, FFT_half, 'r-', linewidth=1.7)
plt.title('FFT of Rect. Pulse', fontsize=21)
plt.xlabel('$f$', fontsize=21, fontweight='bold')
plt.ylabel('S($f$)', fontsize=21, fontweight='bold')
plt.xlim([0,40])
plt.grid()

plt.tight_layout()
plt.show()

Rect.pulse and FFT

As you can see in the picture, the sequence turned out to be acceptable. But the spectrum didn't work out. Increasing the sampling rate doesn't help me at all. What is wrong with my code that causes such an incorrect FFT to be built?

Solution

  • Your discrete Fourier transform is defined as

    enter image description here

    If your f(t) function repeats 3 times in the sequence then you can rewrite this as (neglecting the small difference between N and a multiple of 3!):

    enter image description here

    Now, if m is a multiple of 3 then that last bracketed term is equal to 3, so this Fourier coefficient is (usually) non-zero. If m is not a multiple of 3 then it is 0 (write it out, or sum a geometric series).

    Hence, with 3 repeats you would expect only the 0, 3, 6, 9, 12, 15, … frequency components to be non-zero.

    But, you will say that the 6, 12, 18, … components are zero. However, this is a feature of your particular pulses. In the code below I have adapted the shape of your pulses a bit. Then you get non-zero contributions precisely when m = 0, 3, 6, 9, …

    import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft

A = 3                                                      # amplitude
SR = 1000                                                  # sample size
t = np.linspace( 0, 3, SR, endpoint=False )                # time range from 0 to 3 seconds

n_pulses = 3                                               # number of pulses
tilda = 0.5                                                # pulse duration
t_i = 0.5                                                  # interval between pulses
pulse1_start = 0.25                                        # the beginning of the first pulse

rect_pulses = np.zeros_like( t )
for i in range( n_pulses ):
    pulse_start = pulse1_start + i * ( tilda + t_i )
    pulse_end = pulse_start + tilda
    rect_pulses = np.where( ( t >= pulse_start ) & ( t <= pulse_end ), A - 2.5 * ( t - pulse_start ), rect_pulses )

FFT = fft( rect_pulses )                                   # magnitude of FFT
FFT_n = FFT / len(FFT)                                     # normalization
f_ax = np.linspace( 0, SR, len( FFT ), endpoint=False )    # frequency axis
half_fft = len(FFT) // 2                                   # half the length of the FFT
f_half = f_ax[:half_fft]                                   # frequencies up to f_ax / 2
FFT_half = np.abs(FFT_n[:half_fft])                        # FFT module up to f_ax/2
FFT_half *= 2 # doubling the magnitude 

plt.figure(figsize=(10, 10))
plt.subplot(2, 1, 1)
plt.plot(t, rect_pulses, linewidth=1.7)
plt.title('Pulses', fontsize=21)
plt.xlabel('t', fontsize=21, fontweight='bold')
plt.ylabel('x(t)', fontsize=21, fontweight='bold')
plt.grid()
plt.subplot(2, 1, 2)
plt.plot(f_half, FFT_half, 'r-', linewidth=1.7)
plt.title('FFT of Pulses', fontsize=21)
plt.xlabel('$f$', fontsize=21, fontweight='bold')
plt.ylabel('S($f$)', fontsize=21, fontweight='bold')
plt.xlim([0,40])
plt.grid()

plt.tight_layout()
plt.show()

    enter image description here