简介周期信号的合成与分解实验报告以PDF文件形式提供来源于武汉大学电子信息学院信号与系统类课程实验适合正在学习傅里叶级数、周期信号频谱分析和吉布斯现象的学生阅读。报告内容包括实验目的、基本原理、傅里叶级数性质、所需MATLAB函数说明以及利用有限项级数合成周期对称方波的完整程序与运行结果。程序中分别取前1项、前2项、前5项和前100项进行逼近并另有前5至8项的吉布斯现象观察示例能够帮助理解有限项级数逼近无限项级数时方均误差逐渐减小的过程以及不连续点附近峰起趋于总跳变值约9%的吉布斯现象。资源共1个PDF文件包体大小约608KB图文清晰可直接打印作为实验预习或报告参考。目前已有207人学习下载适合信号与系统、数字信号处理等课程的复习巩固。1. 周期信号傅里叶分解为什么方波合成实验值得用 MATLAB 再跑一遍武汉大学电子信息学院这份《周期信号的合成与分解》实验报告以 PDF 形式在学生群里传了很多年核心内容并不复杂把一个周期方波和三角波拆成傅里叶级数再用 MATLAB 取有限项加回去观察逼近过程和 Gibbs 现象。做信号处理的人对傅里叶级数公式都熟但公式里的系数和有限项截断到底在图形上造成什么影响很多人并没有亲手验证过。这个实验把 50Hz 方波、1/n 和 1/n² 两种收敛速度、FFT 频谱三条线索放在同一组代码里适合作为数字信号处理入门的复现作业也适合工作后用 MATLAB 快速回顾离散谱和吉布斯过冲的本质。下面按实际可运行的顺序拆开讲。2. 方波合成傅里叶级数系数计算与多子图实现2.1 奇对称方波的傅里叶系数推导实验中的方波周期 T0.02s基频 f150Hz角频率 ω100π rad/s。波形关于 t0 奇对称又满足半波对称所以展开式只保留奇次正弦项f(t) Σ (4E/(nπ)) sin(nωt)n1,3,5,...当幅度按实验参数取 E3V总峰峰 6V时基波幅值是 12/π这就是sishu12/pi的来源。需要注意这个变量名容易误导人它并不是每个谐波统一的幅度而是基波系数第 n 次谐波实际是sishu/n。后面所有程序里每一项都写了/(2*i-1)或/(2*i-1)^2就是这个系数的衰减部分。从工程角度看这里有一个值得记住的结论低次谐波系数不会因为后面多取几项而改变。也就是说有限项级数的叠加是单向追加而不是重新分配能量。这也是为什么前 1 项、前 2 项、前 100 项之间可以直接比较而不需要重新归一化。2.2 前 1、2、5、100 项合成程序原始报告里的程序可以整理成下面这样每一段对应 2×2 子图中的一个位置% 奇对称方波合成T0.02sE3V t 0:0.001:0.1; % 0 到 0.1s跨 5 个周期 sishu 12/pi; % 基波系数 4E/pi y sishu*sin(100*pi*t); % 前 1 项只有 50Hz 基波 subplot(221); plot(t, y); axis([0 0.1 -4 4]); xlabel(time); ylabel(前 1 项); y sishu*sin(100*pi*t) sishu*sin(300*pi*t)/3; % 前 2 项 subplot(222); plot(t, y); axis([0 0.1 -4 4]); xlabel(time); ylabel(前 2 项); y sishu*(sin(100*pi*t) sin(300*pi*t)/3 ... sin(500*pi*t)/5 sin(700*pi*t)/7 sin(900*pi*t)/9); subplot(223); plot(t, y); axis([0 0.1 -4 4]); xlabel(time); ylabel(前 5 项); y 0; % 前 100 项用循环累加 for i 1:100 y y sishu*sin((2*i-1)*100*pi*t)/(2*i-1); end subplot(224); plot(t, y); axis([0 0.1 -4 4]); xlabel(time); ylabel(前 100 项);t0:0.001:0.1生成 101 个样点采样间隔 1ms。这个密度画前 5 项时足够平滑但前 100 项的最高谐波是第 199 次频率 9950Hz远高于 1000Hz 采样率对应的 500Hz Nyquist 频率。也就是说高频分量会发生混叠。只是由于这些分量的幅度已经衰减到1/199量级肉眼在 -4V 到 4V 的纵轴上看不出明显异常。如果要做精确仿真先把时间轴改成tlinspace(0, 0.1, 20001)更稳妥。各子图对应关系如下subplot 位置保留项数实际包含的最高谐波最高频率22111 次50 Hz22223 次150 Hz22359 次450 Hz224100199 次9950 Hz这里有一个实验里容易忽略的点前 2 项并不是第 2 次谐波而是第 1、3 次谐波因为方波里偶次项系数为 0。报告中前 N 项指的是 N 个非零项不是前 N 次谐波。2.3 用方均误差做量化收敛验证报告要求观察方均误差随项数增加而减小这个结论可以不用靠感觉直接用 MATLAB 算square_wave 3*sign(sin(100*pi*t)); % 理想方波幅度 ±3V N_list [1 2 5 10 50 100]; for N N_list yN 0; for i 1:N yN yN (12/pi)*sin((2*i-1)*100*pi*t)/(2*i-1); end mse mean((square_wave - yN).^2); fprintf(N%3d, MSE%.6f\n, N, mse); endsign函数把正弦波变成 ±1再乘以 3 得到理想方波。需要说明的是在 t0 处sin(0)0sign(0)0所以理想方波在跳变点有一个 0 值样本这会导致 MSE 略偏大但不会改变随 N 递减的趋势。如果想把跳变点处理得更干净可以用sign(sin(100*pi*(t0.001)))把采样点整体偏移半个步长。这段代码同时验证了另一个结论方均误差不断减小但波形在不连续点处的峰起不会消失这就是下一步 Gibbs 现象要观察的内容。3. Gibbs 现象有限项逼近在间断点处的过冲规律3.1 Gibbs 现象的成因与 9% 过冲傅里叶级数在不连续点附近不是一致收敛的。函数在 t0 处有一个高度为 D 的跳变时部分和的峰起值随项数增加会向跳变点靠拢但峰起高度并不随 N 增大而消失而是趋向约 0.09D。这就是报告里写的约等于总跳变值的 9%。实验中方波从 -3V 跳到 3V跳变值 D6V所以理论上限接近 0.54V。理解这个概念要区分方均误差和最大偏差有限项级数在能量意义下收敛很快但最大偏差在间断点附近收敛极慢。这是两个不同范数下的结论。3.2 前 5、6、7、8 项的程序观察 Gibbs 现象不必重新写代码只把前面合成程序循环化并缩小时间窗到 0.04s两个完整周期让跳变点两侧的细节占据画面% 观察 Gibbs 现象 t 0:0.001:0.04; % 两个周期方便看跳变附近 sishu 12/pi; Ns [5 6 7 8]; for k 1:4 y 0; N Ns(k); for i 1:N y y sishu*sin((2*i-1)*100*pi*t)/(2*i-1); end subplot(2,2,k); plot(t, y); axis([0 0.04 -4 4]); grid on; title(sprintf(前 %d 项, N)); pause(0.3); % 每幅图停 0.3 秒方便对比 endsprintf(前 %d 项, N)用来生成子图标题避免手写五个重复字符串。pause(0.3)在脚本模式下逐幅停留如果想连续播放可以改成pause。四个子图的纵轴统一用axis([0 0.04 -4 4])否则过冲幅度的视觉会被自动缩放掩盖。3.3 过冲怎么读以及常见误区读取过冲不要在整个时间轴上找最大值因为方波平台本身接近 ±3V最大值会落在跳变后的过冲区但需要限定一个窗口jump_t 0.01; % 方波在 0.01s 处从 3 跳到 -3 idx abs(t - jump_t) 0.001; % 取跳变邻域 over max(y(idx)) - 3; % 正方向过冲 fprintf(N%d, 过冲约 %.3f V\n, N, over);这里over max(y(idx)) - 3假设刚好在 3V 平台之后找极值如果跳变方向相反要换成min(y(idx)) 3。窗口大小会影响数值因此报告里一般不写精确到毫伏而是看趋势N峰起位置过冲幅度说明5离跳变点较远约 0.5V 量级波形明显圆滑8更靠近跳变点基本不随 N 下降峰起变窄100几乎贴在跳变点仍接近 0.5V 量级间断处始终存留实际程序里N 从 5 变到 100过冲位置会明显移动但极值大小落在一个很窄的范围内。这个窄范围就是 Gibbs 现象的实验证据。需要留意不要用plot的 LineWidth 或坐标轴缩放去人为放大峰起也不能用smooth或滤波函数把峰起抹掉那已经改变了原始级数。4. 三角波合成与频谱分析1/n² 收敛和 1/n 收敛的差异4.1 偶对称三角波合成程序三角波和方波最大的区别在系数衰减速度。偶对称周期三角波在 T0.02s、幅度按同样峰峰 6V 设置时展开式为直流 3V 加奇次余弦项第 n 项系数按 1/n² 衰减。报告里取的系数sishu24/pi^2配合循环里的/(2*i-1)^2就构成完整谐波幅度。实现代码和方波基本同构% 偶对称周期三角波合成 t 0:0.001:0.1; sishu 24/pi^2; y 3 sishu*cos(100*pi*t); % 前 1 项直流 基波 subplot(221); plot(t, y); axis([0 0.1 -4 4]); xlabel(time); ylabel(前 1 项); y 3 sishu*(cos(100*pi*t) cos(300*pi*t)/9); subplot(222); plot(t, y); axis([0 0.1 -4 4]); xlabel(time); ylabel(前 2 项); y 3 sishu*(cos(100*pi*t) cos(300*pi*t)/9 ... cos(500*pi*t)/25 cos(700*pi*t)/49 ... cos(900*pi*t)/81); subplot(223); plot(t, y); axis([0 0.1 -4 4]); xlabel(time); ylabel(前 5 项); y 3; for i 1:100 n 2*i-1; y y sishu*cos(n*100*pi*t)/n^2; end subplot(224); plot(t, y); axis([0 0.1 -4 4]); xlabel(time); ylabel(前 100 项);n 2*i-1这一步把循环索引换算成实际谐波次数i1 对应 1 次i2 对应 3 次。方波程序里分母是n三角波是n^2这个差异导致三角波前两项就已经能看出大致轮廓而方波前两项仍然是明显的正弦形状。如果波形看起来不对称先检查是不是漏掉了直流分量3。4.2 FFT 频谱分析和频率轴映射报告后半部分用fft画频谱核心代码只有几行N 100; X fft(y, N); % y 取前 100 个采样点 f (1/0.1) * (-N/2 : (N/2-1)); % 频率轴分辨率 10Hz subplot(211); plot(t, y); xlabel(time); ylabel(三角波信号); subplot(212); stem(f, abs(fftshift(X))); xlabel(Frequency(Hz)); ylabel(magnitude);fft(y,N)里 N100 表示取 100 点 DFT(-N/2 : (N/2-1))生成长度 100 的序号向量乘以1/0.110Hz就得到从 -500Hz 到 490Hz 的频率刻度。fftshift把零频移到数组中间舍去负半轴直接用stem(f, abs(X))也可以但频谱形状会被折成左右两半不方便观察。需要提醒的是这里的abs(X)不是真实的信号幅值。MATLAB 的 FFT 结果要还原单边幅值谱通常乘以 2 再除以采样点数 N直流分量则只除以 N。这个实验里只比较方波和三角波各次谐波之间的相对高度统一用abs(fftshift(X))是没问题的但不要拿数值去和12/π直接比。4.3 两信号频谱差异把同一段 FFT 代码分别套在方波和三角波上可以看出明显的谱线衰减差异信号级数形式谐波系数衰减高频成分奇对称方波奇次正弦无直流1/n高频衰减慢丰富偶对称三角波直流奇次余弦1/n²低频占主导高频迅速下降周期信号频谱的三个基本特点在这里都能验证离散性谱线是分立的谐波性每条谱线出现在基频整数倍上收敛性谱线高度随谐波次数增加而递减。方波跳变剧烈所以高频比重大三角波是连续折线所以低频比重大这就是报告里波形变化越剧烈高频成分越多的直接体现。5. 指数形式傅里叶级数复现系数换算与验证技巧5.1 指数形式实现报告思考题要求用指数形式的傅里叶级数重复方波合成。原始程序写成t 0:0.001:0.1; sishu 6/pi; y 0; for i 1:100 n 2*i-1; y y sishu*(exp(1j*n*100*pi*t - 1j*0.5*pi))/n; end plot(t, real(y)); axis([0 0.1 -4 4]); xlabel(time);exp(-j*0.5*pi)等于-j所以这一项相当于-j*6/(πn)*exp(jnωt)。在指数形式的傅里叶级数里这就是正频率方向的系数 F_n。问题在于完整复指数级数要同时累加正负频率正频率一侧如果只加一次再取实部幅值会只有原级数的一半。一种稳妥的改法是y 0; for i 1:100 n 2*i-1; y y sishu*(-1j)*(exp(1j*n*100*pi*t) - exp(-1j*n*100*pi*t))/n; end plot(t, real(y));这里exp(jx)-exp(-jx)2j*sin(x)乘以-j后变成实系数2*sin(x)再乘sishu/n就得到和三角形式一致的结果。如果只保留正频率也可以把系数乘 2但那样不符合指数级数本身的定义也不便于推广到双边谱。5.2 用数值指标验证合成结果做完指数形式复现别只看图用两个指标做闭环ideal 3*sign(sin(100*pi*t)); mse mean((real(y) - ideal).^2); fprintf(MSE %.6f\n, mse); fprintf(最大值 %.3f V\n, max(real(y)));方波的理想幅度是 ±3V所以max(real(y))应该接近3 0.54即过冲后略大于 3V。如果最大值只有 1.5V说明正频率系数差了 2 倍如果 MSE 数量级在 1 以上通常不是 Gibbs 现象而是时间轴采样率或系数漏项的问题。5.3 PDF 报告转代码时的清理习惯这个实验的原始程序经常以 PDF 截图或扫描件形式存在直接从 PDF 复制时减号会变全角j和i会被 OCR 识别成其他字符...续行符在复制后容易丢。常见做法是先在 MATLAB 编辑器里按%注释把程序分成小段逐段运行再检查size(t)和y的长度看到Subscript indices must be real positive integers or logicals这类报错多半是变量i被复数的1j覆盖了。另一个更稳的做法是改动tlinspace(0,0.1,10001)把采样率提到 20kHz 量级再做 FFT这时2*abs(X)/N得到的谐波幅值可以直接和12/π对比验证才真正闭环。本文还有配套的精品资源点击获取
