MATLAB——观察吉布斯效应

吉布斯效应：
当原信号是脉冲信号时，其高频分量主要影响脉冲的跳变沿，低频分量主要影响脉冲的顶部，因此，输出信号波形总是要发生失真，该现象称为吉布斯现象。
原理：
任意周期信号表示为傅立叶级数时，需要无限多项才能逼近原信号，但在实际应用中经常采用有限项级数来代替无限项级数。所选项数越多越接近原信号。所以用一个以原点为中心奇对称的周期性方波，可以用奇次正弦波的叠加来逼近可以观察出吉布斯效应。
代码解释：
主要是下列代码不好理解

N=round((nf+1)/2);
y=zeros(N,max(size(t)));
x=zeros(size(t));
for k=1:2:nf
     x=x+sin(w1*k*t)/k;
     y((k+1)/2,:)=x;
end
mesh(y);

首先弄清变量的含义
nf：要计算的傅立叶级数的项数。
N：正弦函数的采样点数，即离散时间傅立叶变换的长度。
y：一个大小为 $N\times N_t$ 的矩阵，其中 $N_t$ 是时间轴上的采样点数，用于存储正弦函数的傅立叶级数。
x：一个长度为 $N_t$ 的向量，用于存储当前傅立叶级数的计算结果。
t：一个长度为 $N_t$ 的向量，表示时间轴上的采样点。
w1：正弦函数的频率。
k：循环变量，表示当前正在计算的傅立叶级数的项数。
sin(w1kt)：正弦函数的某个傅立叶级数的项。
y((k+1)/2,:)：将当前傅立叶级数的计算结果保存到 y 矩阵中。
mesh(y)：将 y 矩阵绘制成三维图形。
具体来说，代码中的循环从 1 到 nf，每次计算一个正弦函数的傅立叶级数的某一项，并将当前傅立叶级数的计算结果保存到 x 向量中。然后，将 x 向量的前 $(k+1)/2$ 个元素（也就是当前傅立叶级数的计算结果）保存到 y 矩阵中。最后，使用 mesh 函数将 y 矩阵绘制成三维图形。
其中，N=round((nf+1)/2)，这行代码的功能是计算一个整数值 N，它的值是 (nf+1)/2 向最近的整数取整后得到的结果。这个整数 N 的作用是在后面的代码中用来定义一个大小为 N 的向量或矩阵，用于存储计算结果。

具体来说，round 函数是 MATLAB 中的一个内置函数，它的作用是将一个实数值四舍五入为最接近的整数。(nf+1)/2 表示将 (nf+1) 除以 2，计算得到一个实数值，然后使用 round 函数将这个实数值四舍五入为最接近的整数。由于 MATLAB 中的数组索引是从 1 开始的，因此这个整数值 N 的意义是用于定义一个从 1 到 N 的整数序列或向量，用于存储计算结果。

在这段代码中，变量 nf 是用于控制傅立叶级数的项数，因此 N 的值与 nf 的值有关系。N 的值等于 (nf+1)/2 向下取整的结果，这保证了 N 是一个整数。如果 (nf+1)/2 是一个整数，那么 N 的值等于 (nf+1)/2；否则，N 的值比 (nf+1)/2 小 0.5。在这段代码中， (nf+1)/2 的作用是计算一共需要计算多少个傅里叶级数项。因为傅里叶级数中包含正弦和余弦两种函数，所以一共需要计算 (nf+1)/2 个正弦函数和 (nf+1)/2 个余弦函数。如果 nf 是奇数，则正弦函数和余弦函数各需要计算 (nf+1)/2 个；如果 nf 是偶数，则正弦函数和余弦函数各需要计算 nf/2+1 个，其中最后一项为直流分量。但是在实际计算时，我们只需要计算其中的一个，另一个可以通过傅里叶级数的对称性质得到。因此，可以只计算 (nf+1)/2 个正弦函数或 (nf+1)/2 个余弦函数。

y=zeros(N,max(size(t)));
x=zeros(size(t));

这里是生成了两个矩阵。
首先，x=zeros(size(t)) 创建了一个与 t 同样大小的全零矩阵，作为正弦函数的累加器。在后面的循环中，将通过累加一系列正弦函数来生成信号。
其次，y=zeros(N,max(size(t))) 创建了一个大小为 N 行，max(size(t)) 列的全零矩阵，作为输出矩阵。在后面的循环中，将计算 (nf+1)/2 个正弦函数，将它们存储在 y 的每一行中。每一行表示一个正弦函数，而列则表示时间轴上的不同时间点。
因为在循环中需要向 y 中不断添加新的行，所以 y 需要是可变大小的，而 x 则是一个定长的向量，只需要用来累加正弦函数。两个矩阵的大小和类型不同，因此需要分别创建。在这个代码中，我们需要将所有的正弦函数的值存储在 y 矩阵的每一列中，因此需要知道 t 中元素的最大数量。因此，使用 max(size(t)) 可以保证 y 矩阵的列数大于或等于 t 中元素的数量。如果使用 size(t)，则不能保证这一点。单纯的size（t）是指在 MATLAB 中，size(t) 返回一个包含矩阵 t 的行数和列数的向量。如果 t 是一个行向量，则返回值为 [1, N]，其中 N 是 t 中元素的个数。如果 t 是一个列向量，则返回值为 [N, 1]。

for k=1:2:nf

执行一个循环，k 是循环变量，1:2:nf 是循环条件，表示从1开始，以步长为2逐渐递增，直到 nf。

y((k+1)/2,:)=x;

这行代码的含义是将矩阵 x 的值赋值给矩阵 y 中的一行，这一行的索引是 (k+1)/2，而列的索引是 :，表示赋值给整个行。
在这个代码中，循环变量 k 是奇数，因此 (k+1)/2 是一个整数。这个整数表示将 x 中的值赋值给 y 中的第 ((k+1)/2) 行。
因为 x 是一个行向量，所以 y((k+1)/2,:) 将 x 的值复制到 y 的第 ((k+1)/2) 行的所有列中。
总代码：

T1=800;
nf=19;
f1=1/T1;
w1=2*pi*f1;
t=0:1:T1/2;
N=round((nf+1)/2);

y=zeros(N,max(size(t)));
x=zeros(size(t));

for k=1:2:nf;
    x=x+sin(w1*k*t)/k;
    y((k+1)/2,:)=x
end

mesh(y);
axis([0,T1/2,0,N,0,1])

运行结果：