การหาอนุกรมฟูเรียร์ (Fourier Series) ของสัญญาณคาบเวลา และตัวอย่างการเขียนโค้ดด้วย MATLAB และ Python#

---

▷ Fourier Series#

สัญญาณต่อเนื่องเชิงเวลาและเป็นสัญญาณคาบ มีคุณสมบัติดังนี้

$$x(t) = x(t+T), \;\; T > 0, \forall t \in\mathbb{R}$$ โดยที่ $T$ คือ คาบของสัญญาณ (Period)

สัญญาณคาบสามารถเขียนให้อยู่ในรูปแบบของอนุกรมฟูเรียร์ (Fourier Series) ถ้าเป็นไปตามเงื่อนไขของ Dirichlet เช่น ฟังก์ชัน $|x(t)|$ สามารถหาค่าอินทิเกรตได้ในช่วงหนึ่งคาบได้

$$0 < \int_{T} |x(t)| dt < \infty$$

อนุกรมฟูเรียร์ของสัญญาณคาบ $x(t)$ ที่มีความถี่มูลฐาน $f_0$ หรือคาบ $T$ มีรูปแบบดังนี้ (Exponential Complex Form of Fourier Series)

$$x(t) = \sum_{k=-\infty}^{\infty} c_k \, e^{j k\omega_0 t}, \; \omega_0 = 2\pi f_0 = \frac{2\pi}{T}, \; k = 0, \pm 1, \pm 2, ...$$

และสัมประสิทธิ์เชิงซ้อน (Complex Coefficients) ของอนุกรมฟูเรียร์สำหรับสัญญาณ $x(t)$ มีสูตรการคำนวณดังนี้

$$\begin{align} c_0 &= \frac{1}{T} \int_{T} x(t) dt \\ c_k &= \frac{1}{T} \int_{T} x(t)\, e^{-j k\omega_0 t} dt, \; k \neq 0 \end{align}$$

ถ้าใช้สูตรของ Euler

$$e^{\pm j\theta} = cos(\theta) \pm j sin(\theta)$$

สัมประสิทธิ์เชิงซ้อนสามารถเขียนได้ใหม่ ดังนี้ $$c_k = \frac{1}{T} \int_{T} x(t)\, \big[cos(k\omega_0 t) - j sin(k\omega_0 t) \big]\,dt$$

โดยที่

$$\begin{align} \qquad a_k &= \frac{2}{T} \int_{T} x(t)\, cos(k\omega_0 t) dt = a_{-k} \\ b_k &= \frac{2}{T} \int_{T} x(t)\, sin(k\omega_0 t) dt = -b_{-k} \\ c_k &= \frac{1}{2}(a_k - j\, b_k) \\ c^*_k &= \frac{1}{2}(a_k + j\, b_k) = c_{-k} \qquad \mbox{ (Conjugate complex) } \\ a_0 &= c_0 \\ \end{align}$$

ถ้าจะเปลี่ยนจากสัมประสิทธิ์ $c_k$ มาเป็นสัมประสิทธิ์ $a_k$ และ $b_k$ ที่เป็นเลขจำนวนจริง ก็สามารถเขียนฟังก์ชัน $x(t)$ ในรูปแบบของอนุกรมฟูเรียร์ได้ใหม่ดังนี้ (เรียกรูปแบบนี้ว่า Trigonometric Form of Fourier Series)

$$\begin{align} \quad x(t) &= \sum_{k=-\infty}^{\infty} c_k \, e^{j k\omega_0 t} = c_0 + \sum_{k=1}^{\infty} \big( c_{-k}\, e^{-j k\omega_0 t} + c_{k}\, e^{j k\omega_0 t} \big)\\ &= c_0 + \sum_{k=1}^{\infty} \frac{1}{2}(a_k + j b_k)\, e^{-j k\omega_0 t} + \sum_{k=1}^{\infty} \frac{1}{2}(a_k - j b_k)\, e^{j k\omega_0 t} \\ &= a_0 + \sum_{k=1}^{\infty} a_k \frac{(e^{j k\omega_0 t} + e^{-j k\omega_0 t})}{2} - \sum_{k=1}^{\infty} j b_k \frac{(e^{j k\omega_0 t} - e^{-j k\omega_0 t})}{2} \\ &= a_0 + \sum_{k=1}^{\infty} a_k \frac{(e^{j k\omega_0 t} + e^{-j k\omega_0 t})}{2} + \sum_{k=1}^{\infty} b_k \frac{(e^{j k\omega_0 t} - e^{-j k\omega_0 t})}{2j} \\ &= a_0 + \sum_{k=1}^{\infty} a_k\, cos(k\omega_0 t) + \sum_{k=1}^{\infty} b_k\, sin(k\omega_0 t) \\ \end{align}$$

ถ้ากำหนดให้ $$\begin{align} A_0 &= a_0 \\ A_k &= \sqrt{a_k^2 + b_k^2} \\ a_k &= A_k\, cos(\theta_k) \\ b_k &= -A_k\, sin(\theta_k) \\ \theta_k &= -tan^{-1}(\frac{b_k}{a_k}) \end{align}$$

และใช้สูตรตรีโกณมิติต่อไปนี้ $$cos(A+B) = cos(A)cos(B) - sin(A)sin(B)$$

ดังนั้น $x(t)$ จึงเขียนอยู่ในรูปใหม่ดังนี้

$$\begin{align} x(t) &= a_0 + \sum_{k=1}^{\infty} a_k\, cos(k\omega_0 t) + b_k\, sin(k\omega_0 t) \\ &= A_0 + \sum_{k=1}^{\infty} A_k\, cos(\theta_k) cos(k\omega_0 t) - A_k\, sin(\theta_k) sin(k\omega_0 t) \\ x(t) &= A_0 + \sum_{k=1}^{\infty} A_k\, cos(k\omega_0 t + \theta_k) \end{align}$$

 

---

▷ สัญญาณคาบ: Sine & Cosine#

ถ้าเป็นรูปคลื่นไซน์ที่มีคาบเท่ากับ $T$ และ $\omega_0 = 2\pi/T$

$$A\, sin(\omega_0 t) = A \Big(\frac{e^{j\omega_0 t} - e^{-j\omega_0 t} }{2j}\Big)$$

จะได้ค่าสัมประสิทธิ์ $c_k$ (ในกรณีนี้ ไม่จำเป็นต้องใช้การอินทิเกรต)

$$c_k = \begin{cases} +\frac{A}{2j},\; & k = +1 \\ -\frac{A}{2j},\; & k = -1 \\ \;\; 0, \; & |k| \neq 1 \end{cases}$$

ถ้าเป็นฟังก์ชันโคไซน์

$$A\, cos(\omega_0 t) = A \Big( \frac{e^{j\omega_0 t} + e^{-j\omega_0 t} }{2} \Big)$$

จะได้ค่าสัมประสิทธิ์ $c_k$

$$c_k = \begin{cases} \frac{A}{2}, & |k| = 1 \\ 0, & |k| \neq 1 \end{cases}$$

 

---

▷ สัญญาณคาบ: สัญญาณรูปสามเหลี่ยม#

สัญญาณคาบเวลาที่มีลักษณะเป็นรูปสามเหลี่ยม (Triangular Waveform) สามารถเขียนให้อยู่ในรูปแบบของฟังก์ชันทางคณิตศาสตร์ได้ดังนี้

$$x(t) = \begin{cases} \frac{2A}{T}t &, \; 0 \leq t < T/2 \\ 2A(1-\frac{t}{T}) &, \; T/2 \leq t < T \\ x(t \, mod \, T) &, \; \mbox{ otherwise} \end{cases}$$

รูป: สัญญาณรูปสามเหลี่ยม $x(t)$ ในช่วง $t \in [-2T,2T]$

ลองมาดูวิธีการคำนวณหาค่าสัมประสิทธิ์ของอนุกรมฟูเรียร์สำหรับฟังก์ชันนี้

เริ่มต้นด้วยสูตรการอินทิเกรตต่อไปนี้ก่อน

$$\int t\,cos(a\,t)\,dt = \frac{(a\,t)\, sin(a\,t) + cos(a\,t)}{a^2} \\ $$ จากนั้นจึงเริ่มคำนวณ $a_k$ โดยพิจารณาเลือกคาบเวลาอยู่ในช่วง $[-T/2,T/2]$ ซึ่งจะทำได้ง่ายกว่าการคำนวณ $c_k$ สำหรับฟังก์ชันรูปสามเหลี่ยมในกรณีนี้ เนื่องจากเป็นฟังก์ชันคู่

$$\begin{align} a_k &= \frac{2}{T} \int_{-T/2}^{T/2} \frac{2A}{T}|t|\, cos(k\omega_0 t)\, dt \qquad \mbox{ (even function) } \\ &= \frac{4}{T} \int_{0}^{T/2} \frac{2A}{T} t\, cos(k\omega_0 t)\, dt = \frac{8A}{T^2} \int_{0}^{T/2} t\, cos(k\omega_0 t)\, dt \\ &= \frac{8A}{T^2}\Big[ \frac{(k\omega_0\,t)\, sin(k\omega_0\,t) + cos(k\omega_0\,t)}{(k\omega_0)^2} \Big]_{0}^{T/2}\\ &= \frac{8A}{T^2}\Big[ \frac{k\omega_0\,(T/2)\, sin(k\omega_0\,(T/2)) + cos(k\omega_0\,(T/2))}{(k\omega_0)^2} \\ & \qquad \qquad - \frac{k\omega_0\,(0)\, sin(k\omega_0\,(0)) + cos(k\omega_0\,(0))}{(k\omega_0)^2} \Big] \\ &= \frac{8A}{T^2}\Big[ \frac{k\pi\, sin(k\pi) + cos(k\pi)}{(k\omega_0)^2} - \frac{1}{(k\omega_0)^2} \Big] \\ &= \frac{8A}{T^2} \frac{(cos(k\pi)-1)}{(k\omega_0)^2} \\ &= \frac{8A}{T^2} \frac{((-1)^k-1)}{(k\omega_0)^2} \\ &= \frac{8A}{T^2} \frac{((-1)^k-1)}{k^2 (4\pi^2)/T^2} \\ &= 2A\frac{((-1)^k-1)}{k^2 \pi^2} \\ c_k &= \frac{1}{2}(a_k - j b_k) = A\frac{((-1)^k-1)}{k^2 \pi^2} \\ c_0 &= \frac{1}{T} \int_{-T/2}^{T/2} \frac{2A}{T}|t|\,dt = \frac{4A}{T^2} \int_{0}^{T/2} t\,dt = \frac{A}{2} \\ \end{align}$$

และได้ $$x(t) = \frac{A}{2} + \sum_{k=1}^{\infty} A\frac{((-1)^k-1)}{k^2 \pi^2} cos(k\omega_0 t)$$

 

---

▷ ตัวอย่างการเขียนโค้ด MATLAB และสัญญาณรูปสามเหลี่ยม#

ตัวอย่างการเขียนโค้ด MATLAB มีดังนี้ เริ่มต้นด้วยการประกาศใช้ตัวแปรเชิงสัญลักษณ์ (Symbolic Variables) เช่น สัญลักษณ์ $t$ (ตัวแปรอิสระซึ่งเป็นเวลา) $T$ (คาบ) และ $A$ (แอมพลิจูด) จากนั้นสร้างฟังก์ชัน $x_T(t)$ สำหรับหนึ่งคาบของสัญญาณ ซึ่งเป็นฟังก์ชันแบบ piecewise และในตัวอย่างนี้มีลักษณะเป็นรูปสามเหลี่ยม และนำไปใช้ในการสร้างฟังก์ชันที่มีคาบ $x(t)$

เมื่อสร้างฟังก์ชัน $x(t)$ แล้ว จึงนำไปแสดงรูปกราฟของฟังก์ชันโดยใช้คำสั่ง fplot() และแสดงค่าในช่วง $t \in [-2,2]$ และกำหนดค่าตัวเลขให้พารามิเตอร์ $A$ ให้เท่ากับ 1 และ $T$ เท่ากับ 1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clear variables, clear the command window, clear figures
clearvars; clc; clf; close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Code Section 1
syms t real;
syms T A positive;
% Define a triangular-shape function x_T(t) for t in [0,T).
x_T(t) = piecewise( 0 <= t & t < T/2, 2*A*t/T,...
                    T/2 <= t & t < T, 2*A*(1-t/T) );
% Define x(t) as a periodic function.
x(t) = x_T( mod(t,T) );
% Set the specific values for the parameters A and T.
params =  [A,T];
values  = [1,1];
% Set the range for t.
t_range = [-2,2];

% Plot the x(t) function for A=1, T=1 and t in [-2,2]  
figure;
fplot( subs(x(t), params, values), t_range, 'linewidth', 1.5),
axis equal, xlabel('t'), ylabel('x(t)'), grid on;
% Use custom labels.
custom_x_labels = {'-2T', '-T', '0', 'T', '2T'}; 
xticks(linspace( t_range(1), t_range(2), numel(custom_x_labels)));
xticklabels(custom_x_labels);
custom_y_labels = {'0','A'}; 
y_range = [0, eval(subs(A, params, values))];
yticks(linspace( y_range(1), y_range(2), numel(custom_y_labels)))
yticklabels(custom_y_labels);

 

ถัดไปเป็นการคำนวณค่าสัมประสิทธิ์ของอนุกรมฟูเรียร์ $c_k$ รวมถึง $a_k$ และ $b_k$ ด้วย

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Code Section 2
syms k integer;
syms m integer;
omega_0 = 2*pi/T;
% Compute the complex coefficients of the Fourier Series.
c_k = (1/T)*int( x_T(t) * exp(-1j*k*omega_0*t),t,[0,T]);
c_k = simplify( c_k, 'Steps', 100 );
c_k_odd  = simplify( subs(c_k, k, 2*m+1) );
c_k_even = simplify( subs(c_k, k, 2*m ) );
% Compute the a_k coefficients of the Fourier Series.
a_k = simplify( (2/T)*int( x_T(t) * cos(k*omega_0*t),t,[0,T]) );
% Compute the a_0 coefficient of the Fourier Series.
a_0 = (1/T) * int( x_T(t), t, [0,T] );
% Compute the b_k coefficients of the Fourier Series.
b_k = simplify( (2/T)*int( x_T(t) * sin(k*omega_0*t),t,[0,T]) );

% Show the coefficients
c_k, c_k_odd, c_k_even, a_k, b_k, a_0

รูป: ตัวอย่างเอาต์พุตจากการรันโค้ด MATLAB

จากข้อความเอาต์พุต จะได้สัมประสิทธิ์ของอนุกรมฟูเรียร์สำหรับ $x(t)$

$$\begin{align} a_k &= 2A\frac{((-1)^k-1)}{(\pi k)^2}, k > 0 \\ a_k &= -\frac{4A}{(\pi k)^2},\; k=1,3,5,... \\ a_k &= 0,\; k=2,4,6,... \\ a_0 &= c_0 = \frac{A}{2} \\ b_k &= 0 \\ c_k &= \frac{1}{2}(a_k - j b_k) = \frac{a_k}{2} \\ c_k &= -\frac{2A}{(\pi k)^2},\; k=1,3,5,... \\ c_k &= 0,\; k=2,4,6,... \\ \end{align}$$

 

เมื่อได้คำนวณค่าสัมประสิทธิ์ $a_k$ ก็สามารถเขียนฟังก์ชัน $x(t)$ ให้อยู่ในรูปแบบของอนุกรมฟูเรียร์ ถ้าจะลองสร้างฟังก์ชัน $y(t)$ จากอนุกรมฟูเรียร์ของ $x(t)$ สำหรับ $k=0,1,...,K$ ก็มีตัวอย่างดังนี้

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Code Section 3
% Plot a Fourier Series with k=0,1..,K.
K = 5;
y(t) = A/2;
for k=1:1:K
   y(t) = y(t) + 2*A*((-1)^k-1)/(pi*k)^2 * cos(k*omega_0*t);
end
figure;
fplot( subs(y(t),params,values), t_range, 'linewidth', 1.5 ),
axis equal, xlabel('t'), 
ylabel('y(t)'), ylim(y_range), grid on;

รูป: กราฟของฟังก์ชัน $y(t)$

 

ถัดไปเป็นตัวอย่างโค้ดสำหรับวาดรูปกราฟสำหรับค่าสัมประสิทธิ์ $|c_k|$

syms k integer;
c_k = A*((-1)^k-1)/(k*pi)^2; 
K = 10;
ks = 1:1:K;
c_k = [subs(c_k,k,flip(-ks)),A/2,subs(c_k,k,ks)];
figure( 'Position', [100,100,1000,480])
stem( [flip(-ks),0,ks], abs( subs(c_k,A,1) ) ); 
xlabel('k');
ylabel('| c_k |');
title('Magnitude of c_k');
grid on;

รูป: กราฟของฟังก์ชัน $|c_k|$ (Magnitude Plot)

 

---

▷ ตัวอย่างการเขียนโค้ด Python และสัญญาณรูปสามเหลี่ยม#

ถัดไปเป็นตัวอย่างโค้ด Python ใช้ NumPy / SciPy และ SymPy ซึ่งจะให้ผลการทำงานคล้ายกับโค้ดที่เขียนด้วย MATLAB

โค้ดในส่วนแรกเป็นการสร้างฟังก์ชัน $x_T(T)$ และ $x(t)$ ตามลำดับ จากนั้นก็แสดงรูปกราฟของฟังก์ชันดังกล่าว

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt

# Define symbolic variables
t, T, A = sp.symbols('t T A', real=True, positive=True)

# Define a triangular-shape function for t in [0, T)
x_T = sp.Piecewise((2*A*t/T, (0<=t) & (t<T/2)), 
                   (2*A*(1-t/T), (T/2<=t) & (t<T)))

# Define a periodic function
x = x_T.subs(t, t % T)

# Define the t_range
t_range = np.linspace(-2, 2, 400)

fig = plt.figure(figsize=(8, 4))
# Plot the x(t) function for A=1, T=1, and t in [-2,2]
params = {A:1, T:1} 
x_values = [x.subs({t:val}) for val in t_range]
x_values = [x.subs(params) for x in x_values]
plt.plot(t_range, x_values, linewidth=1.5)
plt.axis('equal')
plt.xlabel('t')
plt.ylabel('y(t)')
plt.grid(True)

# Use custom labels
custom_x_labels = ['-2T', '-T', '0', 'T', '2T']
xticks = np.linspace(t_range[0], t_range[-1], len(custom_x_labels))
plt.xticks(xticks, custom_x_labels)
custom_y_labels = ['0', 'A']
y_min = 0
y_max = int(A.subs(params))
yticks = np.linspace(y_min, y_max, len(custom_y_labels))
plt.yticks(yticks, custom_y_labels)
plt.show()

รูป: การวาดรูปกราฟของฟังก์ชัน $x(t)$ ด้วย Python Matplotlib

 

โค้ดส่วนที่สองต่อจากส่วนแรก จะคำนวณสัมประสิทธิ์ของอนุกรมฟูเรียร์สำหรับฟังก์ชัน $x(t)$

# Define symbolic variables
k, t, T = sp.symbols('k t T', integer=True, real=True, positive=True)
omega_0 = 2 * sp.pi / T

x_T = (2*A*t/T*sp.Piecewise((1, t < T/2), (0, True)) 
       + 2*A*(1-t/T)*sp.Piecewise((1, t>= T/2), (0, True)))

# Compute the complex coefficients c_k of the Fourier Series
c_k = sp.simplify( (1/T) * 
      sp.integrate( x_T * sp.exp(-1j*k*omega_0*t), (t,0,T)) )
c_k = sp.simplify( c_k )

# Compute the a_k coefficients of the Fourier Series
a_k = sp.simplify( 
      (2/T) * sp.integrate(x_T*sp.cos(k*omega_0*t), (t,0,T)))

# Compute the a_0 coefficient of the Fourier Series
a_0 = (1/T) * sp.integrate(x_T, (t,0,T))

# Compute the b_k coefficients of the Fourier Series
b_k = sp.simplify( (2/T)*sp.integrate(x_T*sp.sin(k*omega_0*t), (t,0,T)))

# Print the results
print(f"c_k: {c_k}")
print(f"a_k: {a_k}")
print(f"a_0: {a_0}")
print(f"b_k: {b_k}")

ตัวอย่างข้อความเอาต์พุต

c_k: 1.0*A*((-1)**k - 1)/(pi**2*k**2)
a_k: 2*A*((-1)**k - 1)/(pi**2*k**2)
a_0: A/2
b_k: 0

โค้ดส่วนที่สาม สร้างฟังก์ชัน $y(t)$ จากอนุกรมฟูเรียร์ของ $x(t)$ สำหรับ $k=0,1,...,K=11$ แล้วแสดงรูปกราฟที่ได้

# Define parameters
K = 11
A = 1
T = 1
t_range = np.linspace(-2, 2, 400)
omega_0 = 2*np.pi/T
# Initialize y(t) with A/2
y = A/2
# Calculate y(t) for k=1 to K
for k in range(1, K + 1):
    y += 2*A*((-1)**k - 1)/(np.pi*k)**2 * np.cos(k*omega_0*t_range)

# Plot the Fourier Series
plt.plot(t_range, y, linewidth=1.5)
plt.axis('equal')
plt.xlabel('t')
plt.ylabel('y(t)')
plt.grid(True)

# Show the plot
plt.show()

 

---

▷ สัญญาณคาบ: Unit Impulse Train#

ถ้ากำหนดให้สัญญาณ $p(t)$ ได้จากผลรวมอนันต์ของฟังก์ชัน $\delta(t - nT)$ ซึ่งเป็นฟังก์ชัน $\delta(t)$ (มีชื่อว่า Dirac delta impulse) ที่มีการเลื่อนเวลาไปที่ตำแหน่ง $t = nT$ ดังนั้น $p(t)$ จึงเป็นสัญญาณคาบเวลา และมีคาบเท่ากับ $T$

$$p(t) = \sum_{n=-\infty}^{\infty} \delta(t - n T),\; n=0,\pm 1, \pm 2, ...$$

สามารถหาค่าสัมประสิทธิ์ $c_k$ ของอนุกรมฟูเรียร์ ได้ดังนี้

$$\begin{align} c_k &= \frac{1}{T} \int_{T} p(t)\, e^{-j k\omega_0 t} dt \\ &= \frac{1}{T} \int_{-T/2}^{+T/2} \delta(t)\, e^{-j k\omega_0 t} dt = \frac{1}{T} \end{align}$$

ซึ่งจะเห็นได้ว่า $c_k$ ของ $p(t)$ เป็นค่าคงที่และเป็นเลขจำนวนจริง ดังนั้น $p(t)$ จึงเขียนให้อยู่ในรูปแบบของอนุกรมฟูเรียร์ได้ต่อไปนี้ $$\begin{align} p(t) &= \sum_{k=-\infty}^{\infty} c_k \, e^{j k\omega_0 t} = \sum_{k=-\infty}^{\infty} \frac{1}{T} \, e^{j k\omega_0 t}, \;\; \omega_0 = \frac{2\pi}{T} \\ &= \frac{1}{T}\Big[ 1 + \sum_{k=1}^{\infty} e^{-j k \omega_0 t} + e^{j k \omega_0 t} \Big] \\ &= \frac{1}{T}\Big[ 1 + 2 \sum_{k=1}^{\infty} cos(k\omega_0 t) \Big] \\ \end{align}$$

และสรุปได้ว่า

$$p(t) = \sum_{n=-\infty}^{\infty} \delta(t - n T) = \frac{1}{T} \sum_{k=-\infty}^{\infty} e^{j k\omega_0 t} = \frac{1}{T}\Big[ 1 + 2 \sum_{k=1}^{\infty} cos(k\omega_0 t) \Big]$$

คุณสมบัติข้อหนึ่งของ $p(t)$ คือ ผลการอินทิเกรตสำหรับหนึ่งคาบเวลา จะมีค่าเป็น 1

$$\int_{-T/2}^{T/2} p(t)\, dt = \underbrace{\int_{-T/2}^{T/2} \frac{1}{T}\, dt}_{1} + \underbrace{\int_{-T/2}^{T/2} \frac{2}{T} \sum_{k=1}^{\infty} cos(k\omega_0 t)\, dt}_{0} = 1$$

ซึ่งให้ผลเหมือน

$$\int_{-T/2}^{T/2} p(t)\, dt = \int_{-T/2}^{T/2} \delta(t)\, dt = 1$$

 

ถ้าลองเขียนโค้ด MATLAB เพื่อสร้างฟังก์ชัน $x(t)$ ที่มีลักษณะใกล้เคียงกับฟังก์ชัน $p(t)$ โดยใช้ $K$ เป็นพารามิเตอร์ และลองเพิ่มค่าจาก $K=20$ เป็น $K=100$ เพื่อดูแล้วโน้มการเปลี่ยนของรูปกราฟ

$$x(t) = \frac{1}{T}\Big[ 1 + 2 \sum_{k=1}^{K} cos(k\omega_0 t) \Big]$$ $$p(t) = \lim_{K\rightarrow\infty} x(t)$$ ข้อสังเกต เมื่อสร้างฟังก์ชัน $x(t)$ แล้วนำมาแสดงรูปกราฟ จะถูกหารด้วยค่า $K$ ถ้าไม่หารด้วย $K$ เมื่อใช้ค่า $K$ เพิ่มขึ้น แอมพลิจูดของสัญญาณก็จะมีค่ามากขึ้น

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clear variables, clear the command window, clear figures
clearvars; clc; clf; close all;
syms k integer;
syms t real;
syms T positive;
% Set T to 1.
T = subs(T,T,1); 
% Define x(t) as a sum of cosine functions.
x(t) = 1/T;
K=20;
for k=1:1:K
     x(t) = x(t) + (1/T)*cos(k*2*pi/T*t);
end
t_range = [-2,2];
% Plot x(t) by scaling amplitude by 1/K.
figure;
fplot( x(t)/K, t_range), ylim([-0.5,1]), grid on
title(['x(t)/K with K=', num2str(K)])
hold on;
ts = t_range(1):1:t_range(2);
ys = ones(size(ts));
% Create a stem plot
stem(ts, ys, 'r', 'filled'), hold off;

รูป: การวาดกราฟด้วย MATLAB สำหรับฟังก์ชัน $x(t), T=1$ สำหรับ $K=20$ และ $K=100$

 

ถ้าจะเขียนโค้ดด้วย Python + NumPy / SciPy ก็มีตัวอย่างดังนี้

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

# Define symbols and parameters.
t = sp.symbols('t', real=True)

# Set T to 1
T = 1
# Define x(t) as a sum of cosine functions.
x_t = 1/T
# Set the K parameter.
K = 100
omega = 2*sp.pi/T
for k in range(1,K+1):
    x_t += (1/T) * sp.cos(k*omega*t)

# Create a lambda function for x(t) for numerical evaluation.
x = sp.lambdify(t, x_t, 'numpy')
# Define the range for t.
t_range = np.linspace(-2, 2, 1000+1)
plt.figure( figsize=(9,6) )
# Plot x(t) by scaling amplitude by 1/K.
plt.plot(t_range, (1/K)*x(t_range) )
plt.ylim([-0.2, 1.1])
plt.grid(True)
plt.title(f'x(t)/K with K={K}')
# Create a stem plot
ts = np.arange(-2, 3)
ys = np.ones(len(ts))
#plt.stem(ts, ys, linefmt='r', markerfmt='r.', basefmt=' ')
plt.show()

และตัวอย่างรูปกราฟที่ได้

รูป: การวาดกราฟด้วย Python สำหรับฟังก์ชัน $x(t), T=1$ สำหรับ $K=100$

 

---

▷ กล่าวสรุป#

บทความนี้ได้นำเสนอเนื้อหาเกี่ยวกับอนุกรมฟูเรียร์ของสัญญาณคาบเวลา และตัวอย่างการคำนวณหาค่าสัมประสิทธิ์ของอนุกรมฟูเรียร์ โดยการเขียนโค้ดด้วย MATLAB และ Python

 

---

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

Created: 2023-10-15 | Last Updated: 2023-10-18