Fourier Series fit for Triangular pulse in Python

Question

I have a triangular signal (or) more precisely a function that is triangular over a certain interval (0, 20) (I don't care how the output or the input signal is in other ranges). I would like to write the function in terms of cos and sin functions only. Or maybe just trying to find then Fourier series expansion.

So, I have created a function in Python doing this. This, is a symmetric triangular graph over (0, 20) with peak of 400 and minimum of 298.

def Real(t):
    if t<0:
        Ans = 298
    elif t<10:
        Ans = 10.2*t + 298
    elif t<20:
        Ans = 502-10.2*t
    else:
        Ans = 298
    return Ans

ts = np.linspace(0, 20, 100) # for generating plots
Rea = np.array([Real(t) for t in ts])

Then, I got to know about the np.fft.fft(array, n=N) idea (from chatGPT ofcourse). The theory I know is that, basically this command generates $Z_k$ for $ k = 0 \to N $ to satisfy the following equation.

So, I have found the equations and then created another function that would take these coefficients and evaluate the fourier series. And the result was horrible. Below was the code used (for reference)

# Generate time series data over the concerned interval
ts = np.linspace(0, 20, 100)
Rea = np.array([Real(t) for t in ts])

N = 200 # Taking 200 fourier terms
Coeff = np.fft.fft(Rea, n=N)
def Four_Func(t):
    Ans = 0
    for n in range(len(Coeff)):
        Ans += Coeff[n] * np.exp(1j * n * t)
    return np.real(Ans)

Fourier = np.array([Four_Func(t) for t in ts])

The results were plotted as shown

plt.plot(ts, Rea, label='Real Signal in')
plt.plot(ts, Fourier, label='Fourier out')

Answer 1

Your discrete Fourier transform (as in numpy.fft.fft - see its documentation) is defined by

The inverse (accessible as numpy.fft.ifft) is defined by

Hence, if you want to write your series as straight multiples of exponentials

then you will need to scale the original fft; i.e. divide by N:

For a trigonometric series, note that the complex exponentials can be expanded (de Moivre's theorem) as

For a real function, c_m and c_N-m are complex conjugates. The series can be expanded as a pure cosine series if the original function (expanded periodically) is even and a pure sine series if the original function (expanded periodically) is odd. Otherwise, you will need sums of sines and cosines. In this case the periodic extension of your function is even and can be expanded as a cosine series.

You can also get such a series directly from numpy.fft.rfft

A final important point is that periodicity will automatically set f(t=20) equal to f(t=0). You should NOT send this point to the FFT calculator.

Here is code to get both exponential and cosine series for your function. Note that it is only guaranteed to fit at the original points from which you obtained the FFT. For a comment about improving this please see the edits at the bottom of this post.

You can increase NMAX to 100 or higher if you wish, but it will make the graph more congested (and obscure the fact that f(t) should not be explicitly defined at t=20).

import numpy as np
import matplotlib.pyplot as plt

TMAX = 20                                                  # Pulse width
NMAX = 20                                                  # Number of data points (excluding the end, which assumes periodicity)

def RealFunction( t ):                                     # Vectorised version
    return 400 - 10.2 * np.abs( t - 10 )

# Generate time series data over the concerned interval
ts = np.linspace( 0, TMAX, NMAX, endpoint=False )          # IMPORTANT: periodicity will fix the last point - don't set it!
Rea = RealFunction( ts )

N = len( ts )                                              # This is the best you can do with this much data

Coeff = np.fft.fft( Rea )                                  # Formal discrete Fourier transform

# Coefficients for multiples of exponential functions
c = Coeff / N                                              # The 1/N comes strictly from the inverse DFT

# Coefficients for multiples of cosines
b = np.zeros( N//2 )
b[1:] = c[1:N//2] + c[-1:-N//2:-1]                         # Note that exp( ix) = cos(x) + i sin(x)
b[0] = c[0] + c[N//2] * np.cos( np.pi * N / 2 )            #      and  exp(-ix) = cos(x) - i sin(x)

def Exp_Func( t ) :
    Ans = 0
    for m in range( N ):
        Ans += c[m] * np.exp( 2.0 * np.pi * 1j * m * t / TMAX )
    return np.real( Ans )

def Cos_Func( t ) :
    Ans = 0
    for m in range( N//2 ):
        Ans += b[m] * np.cos( 2.0 * np.pi * m * t / TMAX )
    return np.real( Ans )

Exponential_Series = Exp_Func( ts )
Cosine_Series = Cos_Func( ts )

plt.plot( ts, Rea               , 'g-', label="Original" )
plt.plot( ts, Exponential_Series, 'rx', label="Exponential series" )
plt.plot( ts, Cosine_Series     , 'bo', markerfacecolor='none', markeredgecolor='b', label="Cosine series" )
plt.legend()
plt.show()

EDIT: It was raised by OP (KNVCSG) that the fit for the exponential series was not good for intermediate points, whereas the cosine series was much better. The reason for this is that, whereas c[N-m] is the complex conjugate of c[m], exp(2.pi.i.(N-m).t/Tmax) is NOT the complex conjugate of exp(2.pi.i.m.t/Tmax) unless t/Tmax is of the form n/N (i.e., one of the original points).

This is intimately related to the fact that, whilst the DISCRETE Fourier series contains modes going from 0 to N-1, the actual Fourier series in exponentials goes from -INFINITY to +INFINITY. In essence, we would like to "reposition" the second half of the DFT nodes as negative modes, effectively running from -N/2 to +N/2-1. A fudge to accomplish this is to amend the function Exp_func() to the following. When the graph is plotted at points with density 10*NMAX much better agreement is found.

def Exp_Func( t ) :
    Ans = 0
    for m in range( -N // 2, N // 2 ):
        Ans += c[(m+N)%N] * np.exp( 2.0 * np.pi * 1j * m * t / TMAX )
    return np.real( Ans )