ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

傅里叶变换实战:用Python解析方波与三角波频谱

傅里叶变换实战:用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/Nf_s是采样率。很多人画频谱时忘记乘以f_s/N导致横坐标单位错误把“索引”当成“Hz”结果完全对不上理论值。我在第一次用NumPy的fft.fft()函数时就栽过这个跟头——画出来的主瓣位置总在k10怎么算都不对最后发现是忘了做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 EffectFFT只能在f_k k·f_s/N这些离散频率点上“看”频谱。如果真实频率f₀落在两个f_k之间FFT就无法精确捕捉只能显示最近的两个点的能量。提高频率分辨率的唯一办法是增加采样点数N即延长采集时间而不是提高采样率f_s——后者只影响最高可分析频率奈奎斯特频率f_s/2。我曾用一个1024点FFT分析一个100.3 Hz的信号峰出现在k102和k103模糊不清后来改用8192点采集8秒峰就清晰地落在k824附近误差小于0.1 Hz。这说明对于精密频率测量时间就是分辨率。3.3 Python实现的关键参数与物理意义映射用PythonNumPy Matplotlib可视化这些变换核心在于正确建立时域参数与频域参数的映射关系。以下是我反复验证过的标准流程定义时域参数采样率fsHz、信号总时长T秒、采样点数N int(fs * T)。T必须足够长以保证频率分辨率df 1/T足够小例如分析1 Hz信号T至少取1秒。生成时间轴t np.linspace(0, T, N, endpointFalse)。注意endpointFalse避免因端点重复导致的FFT边界效应。构造信号方波square_wave scipy.signal.square(2 * np.pi * f0 * t, duty0.5)。duty0.5确保50%占空比。三角波tri_wave scipy.signal.sawtooth(2 * np.pi * f0 * t, width0.5)。width0.5生成对称三角波。正弦波sin_wave np.sin(2 * np.pi * f0 * t)。计算FFTY np.fft.fft(signal) / N。除以N是关键这使得FFT结果的幅度与信号的时域幅度直接对应Parseval定理。如果不除幅度会随N增大而增大失去物理意义。生成频率轴freqs np.fft.fftfreq(N, 1/fs)。fftfreq()自动处理了正负频率并按[0, f_s/2, ..., -f_s/2df]排序。若只关心正频率取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注意版本兼容性numpy1.21,matplotlib3.5,scipy1.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秒的门函数并计算其连续傅里叶变换的数值近似通过FFTimport 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, endpointFalse) # 以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-, linewidth1.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-, linewidth1.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--, linewidth1.2, label理论 sinc(f*τ)) plt.legend() # 放大主瓣 plt.subplot(2, 2, 3) plt.plot(freqs_rect[pos_mask], np.abs(Y_rect[pos_mask]), r-, linewidth1.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(xfirst_null, colorg, linestyle:, linewidth1.5, labelf理论零点: {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, endpointFalse) f0 5.0 # 基频 5 Hz # 生成方波和三角波 square_wave signal.square(2 * np.pi * f0 * t, duty0.5) tri_wave signal.sawtooth(2 * np.pi * f0 * t, width0.5) # width0.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_collectionTrue, basefmt , markerfmtC0o, linefmtC0-) 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(xh, colorr, linestyle--, alpha0.5) # 三角波频谱 plt.subplot(1, 2, 2) plt.stem(freqs_pos, np.abs(Y_tri_pos), use_line_collectionTrue, basefmt , markerfmtC1o, linefmtC1-) 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(xh, colorr, linestyle--, alpha0.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_collectionTrue是为了提升大数据量下的绘图速度。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, endpointFalse) 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, i1) plt.plot(t, recon, b-, linewidth1.2, labelf{n}项谐波) plt.plot(t, np.sign(np.sin(2*np.pi*f0*t)), r--, linewidth0.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项谐波的重建效果。请聚焦观察每个图的跳变点如t0.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和k01之间。解决方法是使用窗函数如汉宁窗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。例如f05Hz,fs1000Hz则一个周期有200个点。取N200010个周期或N400020个周期就能得到干净的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, endpointFalse) square_wave signal.square(2 * np.pi * f0 * t)这样做之后再看频谱1/n的衰减规律就会清晰呈现。我在调试一个电机控制器时就是因为没遵守这条法则花了两天时间才定位到频谱异常的根源——一个简单的整数倍关系往往就是成败的关键。5.4 “三角波频谱里怎么有偶次谐波理论说只有奇次啊”理论上的理想对称三角波确实只含奇次谐波。但现实中信号发生器的非理想性、ADC的非线性、或代码中的细微偏差都可能引入偶次分量。最常见的原因是scipy.signal.sawtooth()函数在width0.5时生成的是数学上严格的奇函数但如果你在生成t轴时用了endpointTrue会导致首尾点重复破坏了严格的周期性和对称性。终极验证法计算信号的奇偶性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
返回列表