用matlab写一个连续周期信号傅里叶级数

假设我们要求解的连续周期函数为$f(t)$，其周期为$T$，其傅里叶级数为：

$$f(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty} \left(a_n\cos\left(\frac{2\pi n}{T}t\right) + b_n\sin\left(\frac{2\pi n}{T}t\right)\right)$$

其中，

$$a_0 = \frac{2}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} f(t) dt$$

$$a_n = \frac{2}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\cos\left(\frac{2\pi n}{T}t\right) dt$$

$$b_n = \frac{2}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\sin\left(\frac{2\pi n}{T}t\right) dt$$

现在我们可以用 Matlab 实现这个算法。假设我们要求解的函数$f(t)$是一个方波，其周期为$T=2$，其幅值为$1$，在$[-1,0)$上等于$-1$，在$(0,1]$上等于$1$，我们可以这样实现代码：

% 设置参数
T = 2; % 周期
a0 = 0; % a0 的初值
N = 10; % 求和项数

% 计算 a0
a0 = (1/T) * integral(@(t)(-1 + (t>=0)),-T/2,T/2);

% 计算 an 和 bn
an = zeros(N,1);
bn = zeros(N,1);
for n = 1:N
    an(n) = (2/T) * integral(@(t)(-1 + (t>=0)).*cos(2*pi*n*t/T),-T/2,T/2);
    bn(n) = (2/T) * integral(@(t)(-1 + (t>=0)).*sin(2*pi*n*t/T),-T/2,T/2);
end

% 计算傅里叶级数
t = linspace(-T/2,T/2,1000);
f = a0/2;
for n = 1:N
    f = f + an(n)*cos(2*pi*n*t/T) + bn(n)*sin(2*pi*n*t/T);
end

% 绘图
plot(t,f);
xlabel('t');
ylabel('f(t)');
title('傅里叶级数求解方波');

运行上述代码，我们可以得到下图所示的结果：