简介周期信号的合成与分解是信号处理课程中的重要知识点这份PDF文档可作为电子信息类专业学生深入理解傅里叶级数分解与计算机仿真实验的配套参考。文档以武汉大学教学实验报告形式呈现完整覆盖实验目的、基本原理、实验操作与结果分析详细讲解周期信号傅里叶级数的展开式、频谱离散特性、有限项级数逼近方法以及吉布斯现象并给出奇对称方波信号合成所用的MATLAB程序程序包含多种项数取值的对比绘图便于读者对照实验过程理解理论。资源包共1个文件文件类型为PDF压缩包大小约608KB内容精炼且易于保存查阅目前已有207人学习浏览。读者通过该文档可掌握周期信号分解与合成的完整思路了解有限项级数逼近对波形误差的影响同时熟悉频谱分析和实验报告的规范写法适合信号与系统课程实验、期末复习或相关实践项目参考使用。1. 周期信号的合成与分解先把时域和频域的关系摆正一个周期方波在时间轴上是一段一段的跳变但你把它的频谱打出来会发现只有按基频整数倍排列的离散谱线反过来把这些谱线按正确的幅度和相位叠加回去又能把方波重新“长”出来。这就是周期信号的合成与分解最反直觉的地方时域看不连续的信号频域里只是一组数据。它解决的是三个层面的问题用傅里叶级数表达周期信号用 FFT 从采样数据里测出谐波分量以及截断谐波后如何控制振铃和频谱泄漏。适合正在做课程实验、写信号处理前端、或者被“加了窗反而更糊”困扰的工程师。2. 从傅里叶级数到离散频谱分解是投影合成是叠加周期信号分解的数学核心是傅里叶级数但工程上真正会反复用到的形式是复指数级数而不是只含 sin/cos 的三角形式。先把公式和代码对应起来后面用 FFT 时就不会被“为什么要用复数”卡住。2.1 为什么先写复指数级数而不是背正余弦公式设周期为 T、基频为f0 1/T、基波角频率为ω0 2πf0周期信号可以写成x(t) Σ_{k-∞}^{∞} c_k * e^{j k ω0 t}系数通过一个周期内的积分求出来c_k (1/T) ∫_T x(t) * e^{-j k ω0 t} dt这个积分的作用是“投影”因为不同 k 对应的复指数在一个周期内正交积分之后只有 k 那一项留下来其他项贡献为零。正余弦级数需要分别处理 sin、cos 和直流而复指数系数c_k本身就把幅度和相位都装进去了|c_k|是幅度angle(c_k)是初始相位k0对应直流分量。工程上不需要纠结傅里叶级数收敛的充分条件周期信号分段连续、每个周期内能量有限就可以用有限项去逼近。真正要关心的反而是“取多少项”和“采样率够不够”这两点后面会反复出现。2.2 用 Python 对一个周期方波做最小分解用 Python 实现分解时积分直接用矩形法近似采样点间隔dt 1/fs把被积函数在每个采样点上的值乘以dt再累加。下面这段代码生成一个 ±1 的周期方波并求出从-20到20阶的复指数系数。import numpy as np fs 8000 # 采样率 f0 100 # 基频 T 1 / f0 dur 0.5 # 总时长 N int(fs * dur) t np.arange(N) / fs # 前半周期 1后半周期 -1得到周期方波 x np.where((t % T) (T / 2), 1.0, -1.0) # 只取一个完整周期做积分 period_points int(fs * T) ts t[:period_points] xs x[:period_points] K 20 c np.zeros(2 * K 1, dtypecomplex) for idx, k in enumerate(range(-K, K 1)): c[idx] np.sum(xs * np.exp(-2j * np.pi * k * f0 * ts)) * (1 / fs) / T代码里(1 / fs)是矩形积分的dt除以T就是套公式中的1/T。c数组的顺序不是从 1 到 41而是从-20到20所以取第k项时下标要小心。运行后你会发现偶数阶系数接近 0奇数阶系数是纯虚数例如c[21]对应k1约为-0.6366j这正是方波的典型特征。有了系数就能合成回时域。合成是把每个复指数乘以对应系数再相加freqs np.arange(-K, K 1) * f0 x_rec np.zeros_like(ts) for idx, k in enumerate(range(-K, K 1)): x_rec c[idx] * np.exp(2j * np.pi * k * f0 * ts) x_rec np.real(x_rec) print(np.max(np.abs(x_rec - xs)))合成结果与原始方波的最大误差量级在 0.1 到 0.3 之间而且误差集中在跳变点附近。这是因为只用 41 个系数去逼近一个陡峭的方波高频分量被截断了。误差位置比误差大小更有信息量它直接引出后面的吉布斯现象。2.3 幅度谱、相位谱与理论系数对照分解结果不能只看数值最好和理论值对照。周期方波 ±1 的复指数系数理论值为奇数k时c_k -2j / (πk)偶数k时为零。k理论 c_k幅度 |c_k|相位 angle(c_k)00001-2j/π0.6366-π/23-2j/(3π)0.2122-π/25-2j/(5π)0.1273-π/22/4/6000相位清一色是-π/2说明方波主要由正弦分量构成如果信号在时间轴上做了平移相位谱会出现线性旋转但幅度谱不变。这也是周期信号合成与分解的第一个实用结论先看幅度谱决定“有哪些频率”再看相位谱决定“时域波形长什么样”。3. 用 FFT 分解有限长序列分辨率、泄漏和采样对齐理论系数是从连续周期信号积分得到的但实际数据都是采样后的有限长序列。FFT 是最常用的分解工具不过它给出的频率点不是连续的而是按Δf fs/N排列的离散点。这一章解决“测出来的谱和理论谱对不上”的问题。3.1 采样率、点数和频率分辨率的关系FFT 输出第 k 条谱线对应的频率是k * fs / N。相邻两条谱线之间的间隔就是频率分辨率Δf fs / N记录时长N/fs越长分辨率越高。看一个实际参数表采样率 fs点数 N记录时长频率分辨率 Δf800080001s1Hz800040000.5s2Hz1600080000.5s2Hz如果信号里有两个频率只差 0.5Hz就必须让Δf ≤ 0.5Hz否则谱线会连在一起看不出是两个分量。注意提高采样率只提升最高可分析频率不提升分辨率要提升分辨率必须加长记录时长或降低 fs在满足奈奎斯特条件的前提下。3.2 对采样数据做 FFTrfft 与幅度谱的实际写法Python 里对实信号通常用np.fft.rfft它只输出零频到奈奎斯特频率的正频率部分输出长度是N//21。配合np.fft.rfftfreq取频率轴。X np.fft.rfft(x) freq np.fft.rfftfreq(len(x), 1 / fs) # 单边幅度谱正频率谱线乘 2 / N还原时域正弦分量的幅度 amp np.abs(X) * 2 / len(x) amp[0] np.abs(X[0]) / len(x) # 直流分量不乘 2 phase np.angle(X) # 找出幅度最大的几个谱线 top_idx np.argsort(amp)[-5:][::-1] for i in top_idx: print(freq[i], amp[i], phase[i])这里len(x)是本例中的N4000所以Δf 2Hz。方波的 100Hz、300Hz、500Hz 谐波正好落在整数倍频率上每条谱线都对准 bin 中心因此不会泄漏。对amp的缩放要特别注意直流分量只有一个不能乘 2如果N是偶数rfft的最后一个 bin 是奈奎斯特频率同样只有一个也不能乘 2。很多幅度谱看起来偏高或偏低都是栽在这两个 bin 上。3.3 要加窗吗先看是不是整周期采样FFT 隐含对时域序列做周期延拓。如果截取长度不是信号周期的整数倍延拓后会出现人为跳变能量从真实谱线泄漏到相邻频率这就是频谱泄漏。解决办法是加窗但副作用是主瓣变宽。# 故意截一段不完整周期比如 3000 点 bad x[:3000] X_bad np.fft.rfft(bad) # 加汉宁窗后再做 FFT w np.hanning(len(bad)) X_w np.fft.rfft(bad * w) # 幅度恢复时除以窗的累加和不是除以点数 amp_w np.abs(X_w) * 2 / np.sum(w) amp_w[0] np.abs(X_w[0]) / np.sum(w)常用窗的选择可以简单记成矩形窗主瓣最窄但旁瓣高频率分辨能力最好泄漏最明显汉宁窗主瓣宽度约为矩形窗的两倍旁瓣低很多是数值计算默认首选布莱克曼窗旁瓣更低主瓣更宽适合分辨不同频率但不需要很高幅度的场景。窗类型主瓣宽度约第一旁瓣衰减矩形窗2 个 bin-13dB汉宁窗4 个 bin-31dB布莱克曼窗6 个 bin-57dB这些参数的直接后果是如果信号已经满足整周期采样加窗反而会让主瓣变宽、幅值修正变得更复杂如果不满足整周期采样又拒绝加窗就会看到谱线脚下拖出一串平滑下降的“裙边”。所以正确做法是先检查采样时长是不是信号周期的整数倍再决定是否加窗。4. 从方波到吉布斯现象谐波截断会把跳变变成振荡周期方波是最适合观察合成过程的测试信号因为它的频谱最简单但又有无穷多个谐波。有限项合成一定会丢失高频而这个丢失在时域不是简单的“变圆”而是跳变处出现固定幅度的过冲这就是吉布斯现象。4.1 方波系数是 4/πk为什么不能无限逼近跳变点对 ±1 周期方波复指数系数奇数阶是-2j/(πk)。把正负频率合在一起得到实正弦级数x(t) Σ_{k1,3,5,...} (4/πk) * sin(2πk f0 t)每一项的幅度按1/k衰减。按理说取足够多项应该越来越接近方波但方波在跳变点是不连续的傅里叶级数在跳变点的部分和不会均匀收敛而是会在跳变两侧产生幅度不衰减的振荡。这个现象的物理解释是截断高频等价于在频域乘了一个矩形窗时域卷积一个 sinc 函数sinc 函数的侧瓣就是振铃。4.2 峰值不随 K 消失用高频采样把过冲看清楚要观察吉布斯现象采样率必须足够高否则高次谐波会混叠成低频看不到真实的过冲。下面代码把采样率提到 40000Hz这样 101 次谐波的最高频率是 10100Hz仍然低于奈奎斯特频率 20000Hz。def square_partial(K, fs40000, f0100): T 1 / f0 N int(fs * T) t np.arange(N) / fs xr np.zeros(N) # 只累加奇数阶正弦分量 for k in range(1, K 1, 2): xr (4.0 / (np.pi * k)) * np.sin(2 * np.pi * k * f0 * t) return t, xr for K in [5, 21, 101]: t, xr square_partial(K) print(K, xr.max())运行后能看到 K 从 5 增加到 101最大值始终在 1.18 附近比理想上平值 1 高出约 0.18这个过冲量约占跳变幅度 2 的 9%。同时过冲的位置越来越靠近跳变点。这提醒我们单纯靠增加谐波次数并不能消除跳变处的振荡只能在更小的空间尺度上压扁它。4.3 抑制振荡的三种改法改系数而不是只加次数既然矩形窗是振铃的来源就可以用平滑的高频衰减代替硬截断。常见做法是给第 k 个系数乘一个随 k 递减的权重称为 Lanczos sigma 因子fs 40000 f0 100 T 1 / f0 N int(fs * T) t np.arange(N) / fs K 101 x_lanczos np.zeros(N) for k in range(1, K 1, 2): sigma np.sinc(k / (K 1)) # sin(kπ/(K1)) / (kπ/(K1)) x_lanczos (4.0 / (np.pi * k)) * sigma * np.sin(2 * np.pi * k * f0 * t)sigma在 k 接近 K 时衰减到接近 0等效于把频谱边缘慢慢抹掉。相比直接截断过冲会明显降低代价是跳变处的上升沿变缓方波看起来更像梯形。其他的改法包括对系数加汉宁窗形状的包络或者在合成后用低通滤波器处理一次。实际选哪种取决于你要的是最小过冲还是最陡边沿两者无法同时拿到。5. 一健验证合成分解闭环能量守恒和系数回代做完合成与分解不能只看波形像不像还要用数值方法确认信息有没有丢失。Parseval 定理是最快的检查方式。对 DFT 来说时域能量和频域能量满足X_full np.fft.fft(x) energy_time np.sum(np.abs(x) ** 2) energy_freq np.sum(np.abs(X_full) ** 2) / len(x) print(np.abs(energy_time - energy_freq))差值应该在1e-10量级。如果差异很大先检查是不是用了rfft却套用了全 FFT 的归一化公式rfft的半谱不能直接代入这个式子。下一步做系数回代从 FFT 结果里只保留前若干个谐波再逆变换回时域看能量占比是否合理。设delta_f fs / len(x)每个谐波对应的正频率 bin 是harmonic * f0 / delta_fX_full np.fft.fft(x) delta_f fs / len(x) X_keep np.zeros_like(X_full) for harmonic in range(1, 40, 2): # 保留 1,3,5,...,39 次谐波 pos int(round(harmonic * f0 / delta_f)) X_keep[pos] X_full[pos] X_keep[-pos] X_full[-pos] # 负频率部分要同步保留 x_keep np.fft.ifft(X_keep).real ratio np.sum(np.abs(x_keep) ** 2) / np.sum(np.abs(x) ** 2) print(ratio)这里X_keep[pos]对应正频率X_keep[-pos]对应负频率只改一边会导致逆变换结果是复数且幅度减半。ratio 接近 1 说明保留的谐波已经占绝大部分能量如果明显偏低就说明被截断的频率区间里还有不可忽略的能量需要提高保留阶数。自检项预期结果不满足时先查哪里时域频域能量差接近 0FFT 长度、归一化、是否用 rfft合成波形的最大误差随 K 增大而下降采样率是否覆盖谐波范围非整周期采样下的谱线出现泄漏裙边加窗或调整记录时长部分谱回代能量占比随保留谐波数上升负频率 bin 是否漏改本文还有配套的精品资源点击获取 周期信号的合成与分解：从傅里叶级数到FFT实战
简介周期信号的合成与分解是信号处理课程中的重要知识点这份PDF文档可作为电子信息类专业学生深入理解傅里叶级数分解与计算机仿真实验的配套参考。文档以武汉大学教学实验报告形式呈现完整覆盖实验目的、基本原理、实验操作与结果分析详细讲解周期信号傅里叶级数的展开式、频谱离散特性、有限项级数逼近方法以及吉布斯现象并给出奇对称方波信号合成所用的MATLAB程序程序包含多种项数取值的对比绘图便于读者对照实验过程理解理论。资源包共1个文件文件类型为PDF压缩包大小约608KB内容精炼且易于保存查阅目前已有207人学习浏览。读者通过该文档可掌握周期信号分解与合成的完整思路了解有限项级数逼近对波形误差的影响同时熟悉频谱分析和实验报告的规范写法适合信号与系统课程实验、期末复习或相关实践项目参考使用。1. 周期信号的合成与分解先把时域和频域的关系摆正一个周期方波在时间轴上是一段一段的跳变但你把它的频谱打出来会发现只有按基频整数倍排列的离散谱线反过来把这些谱线按正确的幅度和相位叠加回去又能把方波重新“长”出来。这就是周期信号的合成与分解最反直觉的地方时域看不连续的信号频域里只是一组数据。它解决的是三个层面的问题用傅里叶级数表达周期信号用 FFT 从采样数据里测出谐波分量以及截断谐波后如何控制振铃和频谱泄漏。适合正在做课程实验、写信号处理前端、或者被“加了窗反而更糊”困扰的工程师。2. 从傅里叶级数到离散频谱分解是投影合成是叠加周期信号分解的数学核心是傅里叶级数但工程上真正会反复用到的形式是复指数级数而不是只含 sin/cos 的三角形式。先把公式和代码对应起来后面用 FFT 时就不会被“为什么要用复数”卡住。2.1 为什么先写复指数级数而不是背正余弦公式设周期为 T、基频为f0 1/T、基波角频率为ω0 2πf0周期信号可以写成x(t) Σ_{k-∞}^{∞} c_k * e^{j k ω0 t}系数通过一个周期内的积分求出来c_k (1/T) ∫_T x(t) * e^{-j k ω0 t} dt这个积分的作用是“投影”因为不同 k 对应的复指数在一个周期内正交积分之后只有 k 那一项留下来其他项贡献为零。正余弦级数需要分别处理 sin、cos 和直流而复指数系数c_k本身就把幅度和相位都装进去了|c_k|是幅度angle(c_k)是初始相位k0对应直流分量。工程上不需要纠结傅里叶级数收敛的充分条件周期信号分段连续、每个周期内能量有限就可以用有限项去逼近。真正要关心的反而是“取多少项”和“采样率够不够”这两点后面会反复出现。2.2 用 Python 对一个周期方波做最小分解用 Python 实现分解时积分直接用矩形法近似采样点间隔dt 1/fs把被积函数在每个采样点上的值乘以dt再累加。下面这段代码生成一个 ±1 的周期方波并求出从-20到20阶的复指数系数。import numpy as np fs 8000 # 采样率 f0 100 # 基频 T 1 / f0 dur 0.5 # 总时长 N int(fs * dur) t np.arange(N) / fs # 前半周期 1后半周期 -1得到周期方波 x np.where((t % T) (T / 2), 1.0, -1.0) # 只取一个完整周期做积分 period_points int(fs * T) ts t[:period_points] xs x[:period_points] K 20 c np.zeros(2 * K 1, dtypecomplex) for idx, k in enumerate(range(-K, K 1)): c[idx] np.sum(xs * np.exp(-2j * np.pi * k * f0 * ts)) * (1 / fs) / T代码里(1 / fs)是矩形积分的dt除以T就是套公式中的1/T。c数组的顺序不是从 1 到 41而是从-20到20所以取第k项时下标要小心。运行后你会发现偶数阶系数接近 0奇数阶系数是纯虚数例如c[21]对应k1约为-0.6366j这正是方波的典型特征。有了系数就能合成回时域。合成是把每个复指数乘以对应系数再相加freqs np.arange(-K, K 1) * f0 x_rec np.zeros_like(ts) for idx, k in enumerate(range(-K, K 1)): x_rec c[idx] * np.exp(2j * np.pi * k * f0 * ts) x_rec np.real(x_rec) print(np.max(np.abs(x_rec - xs)))合成结果与原始方波的最大误差量级在 0.1 到 0.3 之间而且误差集中在跳变点附近。这是因为只用 41 个系数去逼近一个陡峭的方波高频分量被截断了。误差位置比误差大小更有信息量它直接引出后面的吉布斯现象。2.3 幅度谱、相位谱与理论系数对照分解结果不能只看数值最好和理论值对照。周期方波 ±1 的复指数系数理论值为奇数k时c_k -2j / (πk)偶数k时为零。k理论 c_k幅度 |c_k|相位 angle(c_k)00001-2j/π0.6366-π/23-2j/(3π)0.2122-π/25-2j/(5π)0.1273-π/22/4/6000相位清一色是-π/2说明方波主要由正弦分量构成如果信号在时间轴上做了平移相位谱会出现线性旋转但幅度谱不变。这也是周期信号合成与分解的第一个实用结论先看幅度谱决定“有哪些频率”再看相位谱决定“时域波形长什么样”。3. 用 FFT 分解有限长序列分辨率、泄漏和采样对齐理论系数是从连续周期信号积分得到的但实际数据都是采样后的有限长序列。FFT 是最常用的分解工具不过它给出的频率点不是连续的而是按Δf fs/N排列的离散点。这一章解决“测出来的谱和理论谱对不上”的问题。3.1 采样率、点数和频率分辨率的关系FFT 输出第 k 条谱线对应的频率是k * fs / N。相邻两条谱线之间的间隔就是频率分辨率Δf fs / N记录时长N/fs越长分辨率越高。看一个实际参数表采样率 fs点数 N记录时长频率分辨率 Δf800080001s1Hz800040000.5s2Hz1600080000.5s2Hz如果信号里有两个频率只差 0.5Hz就必须让Δf ≤ 0.5Hz否则谱线会连在一起看不出是两个分量。注意提高采样率只提升最高可分析频率不提升分辨率要提升分辨率必须加长记录时长或降低 fs在满足奈奎斯特条件的前提下。3.2 对采样数据做 FFTrfft 与幅度谱的实际写法Python 里对实信号通常用np.fft.rfft它只输出零频到奈奎斯特频率的正频率部分输出长度是N//21。配合np.fft.rfftfreq取频率轴。X np.fft.rfft(x) freq np.fft.rfftfreq(len(x), 1 / fs) # 单边幅度谱正频率谱线乘 2 / N还原时域正弦分量的幅度 amp np.abs(X) * 2 / len(x) amp[0] np.abs(X[0]) / len(x) # 直流分量不乘 2 phase np.angle(X) # 找出幅度最大的几个谱线 top_idx np.argsort(amp)[-5:][::-1] for i in top_idx: print(freq[i], amp[i], phase[i])这里len(x)是本例中的N4000所以Δf 2Hz。方波的 100Hz、300Hz、500Hz 谐波正好落在整数倍频率上每条谱线都对准 bin 中心因此不会泄漏。对amp的缩放要特别注意直流分量只有一个不能乘 2如果N是偶数rfft的最后一个 bin 是奈奎斯特频率同样只有一个也不能乘 2。很多幅度谱看起来偏高或偏低都是栽在这两个 bin 上。3.3 要加窗吗先看是不是整周期采样FFT 隐含对时域序列做周期延拓。如果截取长度不是信号周期的整数倍延拓后会出现人为跳变能量从真实谱线泄漏到相邻频率这就是频谱泄漏。解决办法是加窗但副作用是主瓣变宽。# 故意截一段不完整周期比如 3000 点 bad x[:3000] X_bad np.fft.rfft(bad) # 加汉宁窗后再做 FFT w np.hanning(len(bad)) X_w np.fft.rfft(bad * w) # 幅度恢复时除以窗的累加和不是除以点数 amp_w np.abs(X_w) * 2 / np.sum(w) amp_w[0] np.abs(X_w[0]) / np.sum(w)常用窗的选择可以简单记成矩形窗主瓣最窄但旁瓣高频率分辨能力最好泄漏最明显汉宁窗主瓣宽度约为矩形窗的两倍旁瓣低很多是数值计算默认首选布莱克曼窗旁瓣更低主瓣更宽适合分辨不同频率但不需要很高幅度的场景。窗类型主瓣宽度约第一旁瓣衰减矩形窗2 个 bin-13dB汉宁窗4 个 bin-31dB布莱克曼窗6 个 bin-57dB这些参数的直接后果是如果信号已经满足整周期采样加窗反而会让主瓣变宽、幅值修正变得更复杂如果不满足整周期采样又拒绝加窗就会看到谱线脚下拖出一串平滑下降的“裙边”。所以正确做法是先检查采样时长是不是信号周期的整数倍再决定是否加窗。4. 从方波到吉布斯现象谐波截断会把跳变变成振荡周期方波是最适合观察合成过程的测试信号因为它的频谱最简单但又有无穷多个谐波。有限项合成一定会丢失高频而这个丢失在时域不是简单的“变圆”而是跳变处出现固定幅度的过冲这就是吉布斯现象。4.1 方波系数是 4/πk为什么不能无限逼近跳变点对 ±1 周期方波复指数系数奇数阶是-2j/(πk)。把正负频率合在一起得到实正弦级数x(t) Σ_{k1,3,5,...} (4/πk) * sin(2πk f0 t)每一项的幅度按1/k衰减。按理说取足够多项应该越来越接近方波但方波在跳变点是不连续的傅里叶级数在跳变点的部分和不会均匀收敛而是会在跳变两侧产生幅度不衰减的振荡。这个现象的物理解释是截断高频等价于在频域乘了一个矩形窗时域卷积一个 sinc 函数sinc 函数的侧瓣就是振铃。4.2 峰值不随 K 消失用高频采样把过冲看清楚要观察吉布斯现象采样率必须足够高否则高次谐波会混叠成低频看不到真实的过冲。下面代码把采样率提到 40000Hz这样 101 次谐波的最高频率是 10100Hz仍然低于奈奎斯特频率 20000Hz。def square_partial(K, fs40000, f0100): T 1 / f0 N int(fs * T) t np.arange(N) / fs xr np.zeros(N) # 只累加奇数阶正弦分量 for k in range(1, K 1, 2): xr (4.0 / (np.pi * k)) * np.sin(2 * np.pi * k * f0 * t) return t, xr for K in [5, 21, 101]: t, xr square_partial(K) print(K, xr.max())运行后能看到 K 从 5 增加到 101最大值始终在 1.18 附近比理想上平值 1 高出约 0.18这个过冲量约占跳变幅度 2 的 9%。同时过冲的位置越来越靠近跳变点。这提醒我们单纯靠增加谐波次数并不能消除跳变处的振荡只能在更小的空间尺度上压扁它。4.3 抑制振荡的三种改法改系数而不是只加次数既然矩形窗是振铃的来源就可以用平滑的高频衰减代替硬截断。常见做法是给第 k 个系数乘一个随 k 递减的权重称为 Lanczos sigma 因子fs 40000 f0 100 T 1 / f0 N int(fs * T) t np.arange(N) / fs K 101 x_lanczos np.zeros(N) for k in range(1, K 1, 2): sigma np.sinc(k / (K 1)) # sin(kπ/(K1)) / (kπ/(K1)) x_lanczos (4.0 / (np.pi * k)) * sigma * np.sin(2 * np.pi * k * f0 * t)sigma在 k 接近 K 时衰减到接近 0等效于把频谱边缘慢慢抹掉。相比直接截断过冲会明显降低代价是跳变处的上升沿变缓方波看起来更像梯形。其他的改法包括对系数加汉宁窗形状的包络或者在合成后用低通滤波器处理一次。实际选哪种取决于你要的是最小过冲还是最陡边沿两者无法同时拿到。5. 一健验证合成分解闭环能量守恒和系数回代做完合成与分解不能只看波形像不像还要用数值方法确认信息有没有丢失。Parseval 定理是最快的检查方式。对 DFT 来说时域能量和频域能量满足X_full np.fft.fft(x) energy_time np.sum(np.abs(x) ** 2) energy_freq np.sum(np.abs(X_full) ** 2) / len(x) print(np.abs(energy_time - energy_freq))差值应该在1e-10量级。如果差异很大先检查是不是用了rfft却套用了全 FFT 的归一化公式rfft的半谱不能直接代入这个式子。下一步做系数回代从 FFT 结果里只保留前若干个谐波再逆变换回时域看能量占比是否合理。设delta_f fs / len(x)每个谐波对应的正频率 bin 是harmonic * f0 / delta_fX_full np.fft.fft(x) delta_f fs / len(x) X_keep np.zeros_like(X_full) for harmonic in range(1, 40, 2): # 保留 1,3,5,...,39 次谐波 pos int(round(harmonic * f0 / delta_f)) X_keep[pos] X_full[pos] X_keep[-pos] X_full[-pos] # 负频率部分要同步保留 x_keep np.fft.ifft(X_keep).real ratio np.sum(np.abs(x_keep) ** 2) / np.sum(np.abs(x) ** 2) print(ratio)这里X_keep[pos]对应正频率X_keep[-pos]对应负频率只改一边会导致逆变换结果是复数且幅度减半。ratio 接近 1 说明保留的谐波已经占绝大部分能量如果明显偏低就说明被截断的频率区间里还有不可忽略的能量需要提高保留阶数。自检项预期结果不满足时先查哪里时域频域能量差接近 0FFT 长度、归一化、是否用 rfft合成波形的最大误差随 K 增大而下降采样率是否覆盖谐波范围非整周期采样下的谱线出现泄漏裙边加窗或调整记录时长部分谱回代能量占比随保留谐波数上升负频率 bin 是否漏改本文还有配套的精品资源点击获取