傅里叶变换实战：用Python解析方波与三角波频谱-育师

1. 这不是数学课，是信号处理的“显微镜”入门

你有没有试过把一段嘈杂的录音里的人声单独提出来？或者在手机拍照时，系统自动抹掉照片里的摩尔纹？又或者，为什么MP3文件比WAV小那么多，听起来却差不多？这些背后，都站着一个叫傅里叶变换的工具——它不是高不可攀的纯数学符号游戏，而是一把实实在在的工程“显微镜”，专门用来拆解任何复杂信号的组成成分。我做嵌入式音频处理项目那会儿，第一次用它把啸叫声从麦克风采集的混合信号里精准切出来，那种“原来声音真的能被拆成一堆正弦波”的震撼感，至今记得清清楚楚。标题里提到的方形函数和三角函数，就是这把显微镜下最典型、也最值得深挖的两个“标本”。方形函数代表的是数字世界里最基础的开关信号——比如单片机IO口输出的高低电平、PWM控制电机的方波、甚至数字通信里的0和1；而三角函数，尤其是正弦和余弦，是构成所有周期性现象的“原子单位”，是傅里叶变换唯一认可的“语言”。这篇笔记不讲ε-δ定义，也不推导积分公式，而是直接带你用Python画出它们的频谱图，亲手验证：一个看起来棱角分明的方波，它的频谱居然由无穷多个奇次谐波叠加而成；一个平滑的三角波，其高频分量衰减得比方波快得多——这种直观对比，才是理解傅里叶变换工程价值的第一步。无论你是电子工程师调试电路噪声，还是算法工程师优化图像滤波器，抑或只是好奇手机语音降噪原理的学生，只要你想搞懂“信号到底由什么组成”，这篇实操笔记就值得你花30分钟跟着敲一遍代码。

2. 为什么选方形函数和三角函数当“试验田”？

2.1 方形函数：数字世界的“硬开关”，频谱藏着无限细节

方形函数（Rectangular Function），在工程中更常被称为“门函数”或“方波”，它的数学定义很简单：在区间[-T/2, T/2]内值为1，其余位置为0。但正是这个看似简单的函数，其傅里叶变换结果却极具教学价值和现实意义。它的频谱是著名的sinc函数：F(ω) = T·sinc(ωT/2)，其中sinc(x) = sin(x)/x。这个表达式背后藏着三个关键事实：第一，它的主瓣宽度与脉冲宽度T成反比——脉冲越窄，频谱越宽，意味着它包含更多高频成分；第二，旁瓣永不消失，且衰减缓慢（按1/ω速度），这解释了为什么数字电路中陡峭的边沿会产生强烈的电磁干扰（EMI）；第三，零点出现在ω = ±2π/T, ±4π/T…处，这些零点位置直接决定了数字信号传输中奈奎斯特采样率的理论下限。我当年调试一款高速ADC采集板时，发现PCB上某段走线像天线一样辐射出强烈干扰，用频谱仪一测，峰值正好落在方波基频的3次、5次谐波上——这正是sinc函数旁瓣的体现。如果当时只关注时域波形的“看起来很干净”，而没去分析它的频域特性，问题根本无从定位。所以，方形函数不是抽象符号，它是数字逻辑电平跳变、开关电源MOSFET导通/关断、甚至雷达脉冲发射的数学化身。它的频谱告诉我们：想让数字信号“干净”，不能只看上升沿够不够快，更要考虑如何抑制那些顽固的高频旁瓣——这直接引出了后续的滤波器设计、PCB布局布线、甚至扩频调制等一整套工程实践。

2.2 三角函数：自然界的“基本音符”，频谱揭示能量分布规律

三角函数在这里特指正弦函数sin(ω₀t)和余弦函数cos(ω₀t)，它们是傅里叶变换的“本征函数”——对它们做傅里叶变换，结果是纯粹的频域冲击函数（Dirac Delta），即能量完全集中在单一频率ω₀上。这个性质看似平凡，却是整个频域分析的基石。现实中不存在完美的单频正弦波，但很多物理现象都高度近似：交流电网电压接近50Hz正弦波，扬声器振膜的微小振动可分解为多个正弦模式，甚至人体心电图（ECG）的R波也能用几个正弦分量拟合。更重要的是，三角函数的线性组合能力，构成了傅里叶级数的核心。一个周期为T的任意函数f(t)，只要满足狄利克雷条件，就能表示为：f(t) = a₀/2 + Σ[aₙcos(nω₀t) + bₙsin(nω₀t)]，其中ω₀ = 2π/T。这里的aₙ和bₙ，就是该函数在n次谐波频率上的“音量旋钮”。我做过一个振动监测项目，用加速度传感器采集电机轴承的振动信号。原始时域波形杂乱无章，但经过FFT（快速傅里叶变换）后，频谱图上清晰出现了基频（对应转速）、2倍频（可能反映不对中）、以及特定的高频峰（指向滚动体缺陷）。这些峰的位置和幅度，直接对应着傅里叶级数中的aₙ、bₙ系数。因此，三角函数不仅是数学工具，更是我们解读物理世界“语言”的字典——每个峰值都在告诉你，系统内部某个部件正在以特定频率、特定强度振动。理解这一点，你就明白为什么机械故障诊断、音频均衡器、甚至MRI成像，都离不开对三角函数基底的运用。

2.3 二者对比：光滑度决定频谱“尾巴”的长短

将方形函数和三角函数放在一起对比，最直观的差异在于它们的光滑性（Smoothness）。方形函数在跳变点处不连续（0阶不连续），其导数在跳变点处是无穷大（1阶不连续）；而三角函数本身及其任意阶导数都是连续且有界的。这个数学性质，在频域上表现为截然不同的能量衰减规律。方形函数的频谱幅度|F(ω)| ~ 1/|ω|，衰减缓慢；而三角函数（作为单频分量）的频谱是δ函数，能量集中；更进一步，一个三角波（由多个正弦波叠加而成的周期函数）的频谱幅度则按1/n²衰减。这意味着，如果你要设计一个低通滤波器来平滑一个方波，让它更接近三角波，你需要的滤波器滚降特性就必须足够陡峭——因为方波的高频能量太多，而三角波的高频能量早已衰减殆尽。我在设计一个LED调光电路时就遇到这个问题：直接用MCU PWM输出方波驱动LED，人眼虽看不出闪烁，但用手机摄像头拍摄会出现明显条纹（摩尔纹），这是因为CMOS传感器的采样频率与PWM频率产生了混叠。后来改用硬件RC滤波器将方波“圆滑”成近似三角波再驱动LED，条纹立刻消失——这本质上就是利用了三角波频谱高频分量更少的特性。所以，选择这两个函数作为学习对象，绝非偶然：它们代表了信号从“最不光滑”到“最光滑”的两个极端，其频谱对比，就是一本活生生的《信号带宽与系统响应》教科书。

3. 核心细节解析：从数学定义到Python可视化

3.1 方形函数的两种形态：非周期门函数 vs 周期方波

在傅里叶分析中，“方形函数”常指两种不同但密切相关的对象：非周期的门函数（Rect(t)）和周期性的方波（Square Wave）。前者是傅里叶变换（FT）的直接对象，后者则是傅里叶级数（FS）的典型示例。门函数Rect(t)定义为：当|t| ≤ 1/2时，Rect(t) = 1；否则为0。它的傅里叶变换是sinc(f) = sin(πf)/(πf)，这是一个连续的频谱。而周期方波，假设周期为T，占空比为50%，其傅里叶级数展开为：Square(t) = (4/π)·[sin(ω₀t) + (1/3)sin(3ω₀t) + (1/5)sin(5ω₀t) + …]，其中ω₀ = 2π/T。这个级数只包含奇次谐波，且幅度按1/n衰减。关键区别在于：FT给出的是连续频谱，适用于非周期信号；FS给出的是离散频谱线，适用于周期信号。实际工程中，我们采集的信号往往是有限长度的，计算机只能处理离散数据，因此必须引入离散傅里叶变换（DFT）及其高效算法FFT。DFT将N个时域采样点x[n]映射到N个频域复数X[k]，其中k代表频率索引。这里有个易错点：DFT的输出X[k]对应的实际频率是f_k = k·f_s/N，f_s是采样率。很多人画频谱时忘记乘以f_s/N，导致横坐标单位错误，把“索引”当成“Hz”，结果完全对不上理论值。我在第一次用NumPy的fft.fft()函数时就栽过这个跟头——画出来的主瓣位置总在k=10，怎么算都不对，最后发现是忘了做f_k = k * fs / N的换算。记住：FFT输出的k值本身没有物理频率意义，必须结合采样率才能还原真实频点。

3.2 三角函数的频谱：δ函数的离散化陷阱

理论上，单频正弦波sin(2πf₀t)的连续傅里叶变换是两个位于±f₀处的δ函数。但在数字世界里，我们永远无法得到真正的δ函数，原因有二：一是有限采样长度（Leakage），二是频率分辨率限制（Resolution）。假设你用采样率f_s = 1000 Hz采集1秒长的sin(2π·100t)信号，共1000个点。理想情况下，FFT应在f = 100 Hz处出现一个尖峰。但若信号频率f₀恰好不是f_s/N的整数倍（例如f₀ = 100.5 Hz），那么能量就会“泄漏”（Spectral Leakage）到相邻的频点上，形成一个展宽的峰，而非尖锐的δ。这是由于DFT隐含地将信号视为周期延拓，当f₀不是基频f_s/N的整数倍时，延拓后会产生不连续，从而引入额外的高频分量。解决方法是使用窗函数（Window Function），如汉宁窗（Hanning），它通过平滑信号两端来减少延拓不连续。另一个陷阱是栅栏效应（Fence Effect）：FFT只能在f_k = k·f_s/N这些离散频率点上“看”频谱。如果真实频率f₀落在两个f_k之间，FFT就无法精确捕捉，只能显示最近的两个点的能量。提高频率分辨率的唯一办法是增加采样点数N（即延长采集时间），而不是提高采样率f_s——后者只影响最高可分析频率（奈奎斯特频率f_s/2）。我曾用一个1024点FFT分析一个100.3 Hz的信号，峰出现在k=102和k=103，模糊不清；后来改用8192点（采集8秒），峰就清晰地落在k=824附近，误差小于0.1 Hz。这说明，对于精密频率测量，时间就是分辨率。

3.3 Python实现的关键参数与物理意义映射

用Python（NumPy + Matplotlib）可视化这些变换，核心在于正确建立时域参数与频域参数的映射关系。以下是我反复验证过的标准流程：

定义时域参数：采样率fs（Hz）、信号总时长T（秒）、采样点数N = int(fs * T)。T必须足够长，以保证频率分辨率df = 1/T足够小（例如分析1 Hz信号，T至少取1秒）。
生成时间轴：t = np.linspace(0, T, N, endpoint=False)。注意endpoint=False，避免因端点重复导致的FFT边界效应。
构造信号：
- 方波：square_wave = scipy.signal.square(2 * np.pi * f0 * t, duty=0.5)。duty=0.5确保50%占空比。
- 三角波：tri_wave = scipy.signal.sawtooth(2 * np.pi * f0 * t, width=0.5)。width=0.5生成对称三角波。
- 正弦波：sin_wave = np.sin(2 * np.pi * f0 * t)。
计算FFT：Y = np.fft.fft(signal) / N。除以N是关键！这使得FFT结果的幅度与信号的时域幅度直接对应（Parseval定理）。如果不除，幅度会随N增大而增大，失去物理意义。
生成频率轴：freqs = np.fft.fftfreq(N, 1/fs)。fftfreq()自动处理了正负频率，并按[0, f_s/2, ..., -f_s/2+df]排序。若只关心正频率，取freqs[:N//2]和Y[:N//2]。
绘制频谱：通常画幅度谱np.abs(Y)。注意，对于实信号，Y是共轭对称的，Y[0]是直流分量，Y[N//2]是奈奎斯特频率分量。

提示：在计算幅度谱前，对Y做np.fft.fftshift()可将零频移到中心，便于观察对称性，但这主要用于教学演示，工程分析中通常用fftfreq的默认顺序。

4. 实操过程：手把手复现频谱图与关键现象

4.1 环境准备与依赖安装

首先，确保你的Python环境已安装核心科学计算库。我推荐使用Anaconda发行版，它预装了大部分所需包。若需手动安装，请执行：

pip install numpy matplotlib scipy

注意版本兼容性：numpy>=1.21,matplotlib>=3.5,scipy>=1.7。旧版本的scipy.signal中某些函数（如sawtooth）行为可能略有不同。我当前使用的是numpy 1.24.3,matplotlib 3.7.1,scipy 1.10.1。创建一个独立的虚拟环境是良好习惯，避免包冲突：

python -m venv fourier_env source fourier_env/bin/activate # Linux/Mac # fourier_env\Scripts\activate # Windows

激活环境后，再运行pip install命令。这样，你的傅里叶实验环境就与系统其他项目完全隔离，结果可复现。

4.2 绘制方形函数（门函数）的傅里叶变换

下面这段代码，将生成一个宽度为2秒的门函数，并计算其连续傅里叶变换的数值近似（通过FFT）：

import numpy as np import matplotlib.pyplot as plt from scipy.fft import fft, fftfreq # 参数设置 fs = 1000 # 采样率 (Hz) T = 4.0 # 总时长 (秒)，需大于门宽，留出零填充空间 N = int(fs * T) # 总采样点数 dt = 1 / fs # 时间步长 # 生成时间轴 t = np.linspace(-T/2, T/2, N, endpoint=False) # 以0为中心，便于观察对称性 # 构造门函数：宽度为2秒，即从-1到+1 rect_width = 2.0 rect_signal = np.where(np.abs(t) <= rect_width/2, 1.0, 0.0) # 计算FFT（注意：FFT假设信号是周期的，所以零填充很重要） # 对rect_signal进行零填充至N点，使其在[-T/2, T/2]内定义 Y_rect = fft(rect_signal) / N # 归一化 freqs_rect = fftfreq(N, dt) # 生成频率轴 # 绘制结果 plt.figure(figsize=(12, 8)) # 时域图 plt.subplot(2, 2, 1) plt.plot(t, rect_signal, 'b-', linewidth=1.5) plt.title('时域：门函数 (Rectangular Function)') plt.xlabel('时间 t (秒)') plt.ylabel('幅度') plt.grid(True) # 频域图（幅度谱） plt.subplot(2, 2, 2) # 只取正频率部分，并转换为Hz pos_mask = freqs_rect >= 0 plt.plot(freqs_rect[pos_mask], np.abs(Y_rect[pos_mask]), 'r-', linewidth=1.5) plt.title('频域：门函数的幅度谱 |F(f)|') plt.xlabel('频率 f (Hz)') plt.ylabel('幅度') plt.grid(True) plt.xlim(0, 5) # 重点关注主瓣区域 # 理论sinc函数对比（红色虚线） f_theory = np.linspace(0, 5, 1000) sinc_theory = np.sinc(f_theory * rect_width) # sinc(x) = sin(πx)/(πx) plt.plot(f_theory, np.abs(sinc_theory), 'k--', linewidth=1.2, label='理论 sinc(f*τ)') plt.legend() # 放大主瓣 plt.subplot(2, 2, 3) plt.plot(freqs_rect[pos_mask], np.abs(Y_rect[pos_mask]), 'r-', linewidth=1.5) plt.title('频域放大：主瓣与第一个零点') plt.xlabel('频率 f (Hz)') plt.ylabel('幅度') plt.grid(True) plt.xlim(0, 1.5) plt.ylim(0, 0.6) # 标注第一个零点 first_null = 1.0 / rect_width # 理论零点在 f = 1/τ = 0.5 Hz plt.axvline(x=first_null, color='g', linestyle=':', linewidth=1.5, label=f'理论零点: {first_null:.1f} Hz') plt.legend() plt.tight_layout() plt.show()

运行这段代码，你会看到四幅图：左上是时域门函数，右上是其频谱（红色实线）与理论sinc函数（黑色虚线）的对比，左下是主瓣放大图，并用绿色虚线标出了理论第一个零点位置（f = 1/τ = 0.5 Hz）。关键观察点：实线与虚线几乎完全重合，证明了数值FFT对连续FT的良好逼近；主瓣宽度（从第一个零点到第二个零点）约为1 Hz，这与理论Δf ≈ 1/τ完美吻合。这个实验直接验证了“时域越窄，频域越宽”的核心原理。

4.3 绘制周期方波与三角波的傅里叶级数对比

接下来，我们对比周期方波和三角波的频谱，重点观察谐波衰减规律：

import numpy as np import matplotlib.pyplot as plt from scipy import signal # 参数设置 fs = 1000 T = 2.0 # 采集2秒，包含多个完整周期 N = int(fs * T) t = np.linspace(0, T, N, endpoint=False) f0 = 5.0 # 基频 5 Hz # 生成方波和三角波 square_wave = signal.square(2 * np.pi * f0 * t, duty=0.5) tri_wave = signal.sawtooth(2 * np.pi * f0 * t, width=0.5) # width=0.5 生成对称三角波 # 计算FFT Y_square = np.fft.fft(square_wave) / N Y_tri = np.fft.fft(tri_wave) / N freqs = np.fft.fftfreq(N, 1/fs) # 只取正频率 pos_mask = freqs >= 0 freqs_pos = freqs[pos_mask] Y_square_pos = Y_square[pos_mask] Y_tri_pos = Y_tri[pos_mask] # 绘制 plt.figure(figsize=(14, 6)) # 方波频谱 plt.subplot(1, 2, 1) plt.stem(freqs_pos, np.abs(Y_square_pos), use_line_collection=True, basefmt=" ", markerfmt="C0o", linefmt="C0-") plt.title('方波的幅度谱') plt.xlabel('频率 (Hz)') plt.ylabel('幅度') plt.grid(True) plt.xlim(0, 50) plt.ylim(0, 0.5) # 标注奇次谐波 harmonics_square = [f0 * n for n in range(1, 11, 2)] # 1,3,5,...,19次 for h in harmonics_square: if h <= 50: plt.axvline(x=h, color='r', linestyle='--', alpha=0.5) # 三角波频谱 plt.subplot(1, 2, 2) plt.stem(freqs_pos, np.abs(Y_tri_pos), use_line_collection=True, basefmt=" ", markerfmt="C1o", linefmt="C1-") plt.title('三角波的幅度谱') plt.xlabel('频率 (Hz)') plt.ylabel('幅度') plt.grid(True) plt.xlim(0, 50) plt.ylim(0, 0.15) # 标注奇次谐波（三角波也是奇次谐波为主） harmonics_tri = [f0 * n for n in range(1, 11, 2)] for h in harmonics_tri: if h <= 50: plt.axvline(x=h, color='r', linestyle='--', alpha=0.5) plt.tight_layout() plt.show()

运行此代码，你会得到两张并排的频谱图。关键观察与分析：

方波图（左）：清晰可见1、3、5、7、9...次谐波的尖峰，幅度依次为1、1/3、1/5、1/7、1/9...，完美符合1/n衰减律。在50Hz范围内，你能看到9次谐波（45Hz），其幅度约为基频的1/9≈0.11，与图中高度一致。
三角波图（右）：同样有1、3、5、7、9...次谐波，但幅度衰减快得多。理论值应为1/n²，即1、1/9、1/25、1/49、1/81...。图中9次谐波（45Hz）的幅度已非常微弱（<0.015），远低于方波同次谐波。这直观解释了为何三角波听起来比方波“柔和”——它的高频能量被大幅抑制。

注意：stem()绘图比plot()更能体现离散频谱的特性。use_line_collection=True是为了提升大数据量下的绘图速度。

4.4 揭示“吉布斯现象”：方波重建中的永恒振铃

傅里叶级数的一个著名现象是吉布斯现象（Gibbs Phenomenon）：当用有限项正弦波叠加来逼近一个有跳变的函数（如方波）时，在跳变点附近会出现一个不随项数增加而消失的过冲（Overshoot），其峰值约为跳变幅度的9%。这是由sinc函数的旁瓣能量造成的，是数学上的必然，而非计算误差。下面代码演示这一现象：

import numpy as np import matplotlib.pyplot as plt def reconstruct_square_wave(t, f0, n_terms): """用n_terms个奇次谐波重建方波""" omega0 = 2 * np.pi * f0 wave = np.zeros_like(t) for n in range(1, n_terms*2, 2): # 取1,3,5,...,2n_terms-1 wave += (4 / (np.pi * n)) * np.sin(n * omega0 * t) return wave # 参数 f0 = 1.0 t = np.linspace(0, 4, 2000, endpoint=False) n_list = [1, 5, 20, 100] plt.figure(figsize=(12, 8)) for i, n in enumerate(n_list): recon = reconstruct_square_wave(t, f0, n) plt.subplot(2, 2, i+1) plt.plot(t, recon, 'b-', linewidth=1.2, label=f'{n}项谐波') plt.plot(t, np.sign(np.sin(2*np.pi*f0*t)), 'r--', linewidth=0.8, label='理想方波') # 理想方波参考 plt.title(f'吉布斯现象：{n}项傅里叶级数重建') plt.xlabel('时间 t (秒)') plt.ylabel('幅度') plt.grid(True) plt.ylim(-1.2, 1.2) plt.legend() plt.tight_layout() plt.show()

运行后，四张子图展示了从1项到100项谐波的重建效果。请聚焦观察每个图的跳变点（如t=0.5, 1.0, 1.5...）：即使到了100项，跳变点处依然存在约0.09的过冲和下冲，且这个过冲的宽度随着项数增加而变窄，但高度恒定。这深刻说明：任何试图用有限带宽信号完美复现跳变的努力，都会在跳变处留下一个无法消除的“振铃”。在工程中，这提醒我们：视频编码中的块效应、音频压缩中的预回声（Pre-echo），其根源都与此相关。解决方案不是无限增加带宽，而是采用更平滑的窗函数或改进的编码算法来抑制旁瓣。

5. 常见问题与排查技巧实录

5.1 “我的频谱图怎么是平的？一点峰都没有！”

这是新手最常见的报错，90%的原因是信号幅度太小或直流偏移过大。FFT对绝对幅度敏感，如果信号本身是微伏级的传感器输出，而你的ADC参考电压是3.3V，那么量化后的数字值可能全为0或只有几个LSB在跳动，FFT自然看不到有效频谱。反之，如果信号有一个很大的直流分量（比如运放输出的2.5V偏置），它会占据FFT结果的Y[0]（直流分量），而其他交流分量的幅度相比之下微不足道，导致频谱图看起来像一条直线。

排查步骤：

先看时域：用plt.plot(t, signal)确认信号在时域是否“动起来”了。如果是一条水平线，检查信号源、接线、ADC配置。
检查直流分量：计算np.mean(signal)。如果远不为0（例如>0.1），说明存在显著直流偏移。解决方法：在FFT前做signal_centered = signal - np.mean(signal)，即“去直流”。
检查幅度范围：print(np.min(signal), np.max(signal))。如果范围极小（如-0.001到0.001），说明信号太弱。需要前置放大或调整ADC增益。
验证FFT输出：打印np.abs(Y)[0]（直流）和np.abs(Y)[10]（假设10Hz处有信号）的值。如果两者相差1000倍以上，说明交流分量被直流淹没。

实操心得：在嵌入式系统中，我习惯在采集后立即做两件事：1)signal -= np.mean(signal)；2)signal /= np.std(signal)（归一化）。这能确保FFT结果具有可比性，避免因量纲不同导致的误判。

5.2 “频谱峰的位置对不上理论值，差了一倍！”

这几乎总是采样率（fs）或频率轴计算错误导致的。最常见的错误是混淆了fftfreq的返回值。np.fft.fftfreq(N, d)中的d是采样间隔（秒），即1/fs，而不是fs本身。写成fftfreq(N, fs)是致命错误，会导致频率轴被放大fs²倍。

快速验证法：

理论基频为f0的正弦波，其FFT峰值应出现在索引k0 = round(f0 * N / fs)处。
手动计算：k0 = int(f0 * N / fs)，然后检查np.abs(Y[k0])是否为最大值之一。
如果k0计算结果是200，但峰值在400，那一定是d参数写反了（用了fs而非1/fs）。

另一个原因是信号频率不是基频的整数倍，导致能量泄漏。此时，峰值不会精确落在k0，而是在k0和k0+1之间。解决方法是使用窗函数（如汉宁窗）：windowed_signal = signal * np.hanning(N)，然后再做FFT。窗函数会牺牲一点频率分辨率（主瓣变宽），但能极大抑制旁瓣，使主峰更清晰。

5.3 “为什么我的方波频谱没有‘1/n’衰减？看起来像随机噪声！”

这通常是采样点数N不够，或信号周期未被整除造成的。FFT要求信号在N点内是“周期性”的。如果方波的周期T0与总时长T不匹配（即T/T0不是整数），那么FFT会将信号视为在T处突然截断，产生不连续，从而引入大量虚假的高频分量，掩盖了真实的谐波结构。

黄金法则：确保N是方波周期采样点数的整数倍。即，设方波周期T0 = 1/f0，则N应满足N % (f_s * T0) == 0。例如，f0=5Hz,fs=1000Hz，则一个周期有200个点。取N=2000（10个周期）或N=4000（20个周期），就能得到干净的1/n谱。

快速修复代码：

# 计算一个周期的采样点数 points_per_cycle = int(fs / f0) # 调整N为points_per_cycle的整数倍 N = points_per_cycle * 10 # 取10个完整周期 t = np.linspace(0, N/fs, N, endpoint=False) square_wave = signal.square(2 * np.pi * f0 * t)

这样做之后，再看频谱，1/n的衰减规律就会清晰呈现。我在调试一个电机控制器时，就是因为没遵守这条法则，花了两天时间才定位到频谱异常的根源——一个简单的整数倍关系，往往就是成败的关键。

5.4 “三角波频谱里怎么有偶次谐波？理论说只有奇次啊！”

理论上的理想对称三角波确实只含奇次谐波。但现实中，信号发生器的非理想性、ADC的非线性、或代码中的细微偏差，都可能引入偶次分量。最常见的原因是：scipy.signal.sawtooth()函数在width=0.5时，生成的是数学上严格的奇函数，但如果你在生成t轴时用了endpoint=True，会导致首尾点重复，破坏了严格的周期性和对称性。

终极验证法：

计算信号的奇偶性：odd_part = (signal - signal[::-1]) / 2。理想奇函数的odd_part应等于原信号，even_part = (signal + signal[::-1]) / 2应接近零。
如果np.max(np.abs(even_part)) > 1e-10，说明存在偶分量。此时，强制修正：signal_odd = (signal - np.flip(signal)) / 2，再对此signal_odd做FFT。

此外，检查你的f0和fs是否精确。浮点数计算中，1000/5可能不是精确的200，而是199.999999，导致周期不严格闭合。使用np.arange(N) * dt代替np.linspace可以避免此类问题。

实操心得：在做精密频谱分析时，我从不依赖linspace生成时间轴。一律用t = np.arange(N) * (1/fs)，因为它保证了dt的绝对精度，消除了浮点累积误差。这个小习惯，帮我避开了无数个“幽灵谐波”的坑。

6. 工程延伸：从笔记到真实世界的信号处理

6.1 从方波频谱到EMI滤波器设计

理解方波的sinc频谱，是设计电源输入滤波器的第一步。一个5V/1A的DC-DC降压芯片，其开关节点（SW）的方波电压，基频可能是500kHz。根据sinc函数，其显著能量会延伸到1/τ之外，其中τ是开关管的上升/下降时间。若τ = 10ns，则1/τ = 100MHz——这意味着，你的滤波器必须在100MHz以内都有效，否则高频噪声会通过空间辐射出去，导致产品EMC测试失败。实践中，我们会用LC滤波器，其截止频率f_c = 1/(2π√(LC))。但仅设f_c = 1MHz是远远不够的，因为sinc的旁瓣在10MHz、100MHz仍有能量。因此，工程师会采用多级滤波：一级LC抑制基频和谐波，二级π型滤波器（L-C-L）进一步衰减高频，再配合铁氧体磁珠吸收超高频。所有这些决策，源头都来自对方波频谱的深刻理解。我曾参与一个医疗设备项目，初版设计EMC