数学建模优秀论文算法-傅里叶变化

傅里叶变换小白教程：从“拆解音乐”到“解读世界”

1. 引言：复杂事物的“简单密码”

你有没有想过：

• 一首钢琴曲为什么能被拆解成“do、re、mi”等简单音符？
• 一幅油画为什么能通过“红、绿、蓝”三原色混合而成？
• 一段电流信号为什么能去掉噪声、保留有用信息？

答案藏在傅里叶变换里——它是一把“拆解复杂事物的钥匙”，能将任何满足条件的信号（声音、图像、电流等）分解为不同频率的正弦波/余弦波的叠加，帮我们从“时间的维度”跳转到“频率的维度”，看清事物的本质。

2. 背景溯源：从“热传导”到“傅里叶级数”

傅里叶变换的起源，要从19世纪法国数学家让-巴普蒂斯·约瑟夫·傅里叶（Jean-Baptiste Joseph Fourier）的研究说起。

1807年，傅里叶在解决热传导问题时发现：任何周期函数都可以分解为正弦波和余弦波的叠加。这一结论在当时引发了争议（拉格朗日等数学家认为“非光滑函数无法用光滑的正弦波叠加”），但傅里叶通过严格推导证明了其正确性。1822年，他在《热的解析理论》一书中系统提出傅里叶级数（Fourier Series），为后来的傅里叶变换奠定了基础。

20世纪后，随着通信、信号处理等领域的发展，傅里叶级数被扩展到非周期函数，形成了我们今天熟知的傅里叶变换（Fourier Transform）。

3. 核心思想：所有信号都是“正弦波的叠加”

傅里叶变换的核心可以用一句话概括：

任何满足条件的信号，都能分解为不同频率的正弦波（或复指数波）的线性叠加；反之，这些正弦波也能重构出原信号。

举个最直观的例子：我们听到的音乐是时域信号（随时间变化的空气振动），而音乐的“音调”对应频率（do=261.6Hz，re=293.7Hz，mi=329.6Hz……）。傅里叶变换做的事，就是把“时域的曲子”转换成“频域的乐谱”——告诉你：

• 每个频率的声音有多强（振幅谱）；
• 每个频率的声音在什么时候出现（相位谱）。

4. 基础铺垫：理解傅里叶的“语言”

在正式学习傅里叶变换前，需要先掌握两个关键概念——周期函数和正交函数。

4.1 周期函数与角频率

如果一个函数满足f(t+T)=f(t)f(t + T) = f(t)f(t+T)=f(t)（TTT是常数，称为周期），则它是周期函数。例如：

• 正弦波sin⁡(ω0t)\sin(\omega_0 t)sin(ω0​t)的周期T=2πω0T = \frac{2\pi}{\omega_0}T=ω0​2π​；
• 其中ω0=2πT\omega_0 = \frac{2\pi}{T}ω0​=T2π​称为角频率（单位：rad/s），对应“单位时间内的角度变化”；
• 频率f0=1T=ω02πf_0 = \frac{1}{T} = \frac{\omega_0}{2\pi}f0​=T1​=2πω0​​（单位：Hz，赫兹），表示“每秒振动的次数”。

4.2 正交函数：分解的“坐标系”

要分解一个函数，需要一组正交函数作为“基”——就像用x轴、y轴描述平面上的点。

两个函数g(t)g(t)g(t)和h(t)h(t)h(t)在区间[a,b][a, b][a,b]上正交，当且仅当：∫abg(t)h(t)dt=0(geqh)\int_a^b g(t) h(t) dt = 0 \quad (geq h)∫ab​g(t)h(t)dt=0(geqh)而当g=hg = hg=h时，积分结果是它们的“模长平方”（类似坐标系轴的长度）。

对于周期为TTT的函数，最常用的正交基是三角函数系：{1,cos⁡(ω0t),sin⁡(ω0t),cos⁡(2ω0t),sin⁡(2ω0t),… }\{1, \cos(\omega_0 t), \sin(\omega_0 t), \cos(2\omega_0 t), \sin(2\omega_0 t), \dots\}{1,cos(ω0​t),sin(ω0​t),cos(2ω0​t),sin(2ω0​t),…}这个函数系在[−T/2,T/2][-T/2, T/2][−T/2,T/2]上满足正交性，比如：

• ∫−T/2T/2cos⁡(mω0t)sin⁡(nω0t)dt=0\int_{-T/2}^{T/2} \cos(m\omega_0 t) \sin(n\omega_0 t) dt = 0∫−T/2T/2​cos(mω0​t)sin(nω0​t)dt=0（余弦与正弦正交）；
• ∫−T/2T/2cos⁡(mω0t)cos⁡(nω0t)dt=0\int_{-T/2}^{T/2} \cos(m\omega_0 t) \cos(n\omega_0 t) dt = 0∫−T/2T/2​cos(mω0​t)cos(nω0​t)dt=0（meqnm eq nmeqn时，余弦项之间正交）。

5. 傅里叶级数：周期信号的“拆解工具”

傅里叶级数是傅里叶变换的“前身”，用于处理周期信号（如方波、正弦波）。

5.1 三角形式：直观的“正弦波叠加”

对于周期为TTT、角频率为ω0=2πT\omega_0 = \frac{2\pi}{T}ω0​=T2π​的函数f(t)f(t)f(t)，若满足狄利克雷条件（后面会讲），则可以展开为：f(t)=a02+∑n=1∞[ancos⁡(nω0t)+bnsin⁡(nω0t)]f(t) = \frac{a_0}{2} + \sum_{n=1}^\infty \left[ a_n \cos(n\omega_0 t) + b_n \sin(n\omega_0 t) \right]f(t)=2a0​​+∑n=1∞​[an​cos(nω0​t)+bn​sin(nω0​t)]其中：

• a02\frac{a_0}{2}2a0​​：直流分量（常数项，对应频率0的成分）；
• ancos⁡(nω0t)a_n \cos(n\omega_0 t)an​cos(nω0​t)：余弦谐波，bnsin⁡(nω0t)b_n \sin(n\omega_0 t)bn​sin(nω0​t)：正弦谐波；
• n=1,2,…n = 1, 2, \dotsn=1,2,…：谐波次数——n=1n=1n=1是基波（频率与原函数相同），n>1n>1n>1是高次谐波（频率为基波的整数倍）。

5.2 系数计算：正交性的“魔法”

要得到a0,an,bna_0, a_n, b_na0​,an​,bn​，需利用三角函数系的正交性——让其他项在积分时“消失”，只保留目标系数。

（1）计算直流分量a0a_0a0​

将傅里叶级数两边在[−T/2,T/2][-T/2, T/2][−T/2,T/2]上积分：∫−T/2T/2f(t)dt=∫−T/2T/2a02dt+∑n=1∞[an∫−T/2T/2cos⁡(nω0t)dt+bn∫−T/2T/2sin⁡(nω0t)dt]\int_{-T/2}^{T/2} f(t) dt = \int_{-T/2}^{T/2} \frac{a_0}{2} dt + \sum_{n=1}^\infty \left[ a_n \int_{-T/2}^{T/2} \cos(n\omega_0 t) dt + b_n \int_{-T/2}^{T/2} \sin(n\omega_0 t) dt \right]∫−T/2T/2​f(t)dt=∫−T/2T/2​2a0​​dt+∑n=1∞​[an​∫−T/2T/2​cos(nω0​t)dt+bn​∫−T/2T/2​sin(nω0​t)dt]由于余弦和正弦在整周期内的积分是0，右边求和项全为0，因此：a0=2T∫−T/2T/2f(t)dta_0 = \frac{2}{T} \int_{-T/2}^{T/2} f(t) dta0​=T2​∫−T/2T/2​f(t)dt

（2）计算余弦系数ana_nan​

将傅里叶级数两边乘以cos⁡(nω0t)\cos(n\omega_0 t)cos(nω0​t)并积分：∫−T/2T/2f(t)cos⁡(nω0t)dt=an⋅T2\int_{-T/2}^{T/2} f(t) \cos(n\omega_0 t) dt = a_n \cdot \frac{T}{2}∫−T/2T/2​f(t)cos(nω0​t)dt=an​⋅2T​（其他项因正交性消失）因此：an=2T∫−T/2T/2f(t)cos⁡(nω0t)dt(n≥1)a_n = \frac{2}{T} \int_{-T/2}^{T/2} f(t) \cos(n\omega_0 t) dt \quad (n \geq 1)an​=T2​∫−T/2T/2​f(t)cos(nω0​t)dt(n≥1)

（3）计算正弦系数bnb_nbn​

类似地，乘以sin⁡(nω0t)\sin(n\omega_0 t)sin(nω0​t)并积分：bn=2T∫−T/2T/2f(t)sin⁡(nω0t)dt(n≥1)b_n = \frac{2}{T} \int_{-T/2}^{T/2} f(t) \sin(n\omega_0 t) dt \quad (n \geq 1)bn​=T2​∫−T/2T/2​f(t)sin(nω0​t)dt(n≥1)

5.3 例子：方波的傅里叶级数分解

我们用周期方波验证傅里叶级数的效果。假设方波周期T=2T = 2T=2（ω0=π\omega_0 = \piω0​=π），表达式为：f(t)={1∣t∣<0.500.5<∣t∣<1f(t) = \begin{cases} 1 & |t| < 0.5 \\ 0 & 0.5 < |t| < 1 \end{cases}f(t)={10​∣t∣<0.50.5<∣t∣<1​

计算系数：

• a0=∫−11f(t)dt=1a_0 = \int_{-1}^{1} f(t) dt = 1a0​=∫−11​f(t)dt=1（直流分量）；
• an=∫−0.50.5cos⁡(nπt)dt=2nπsin⁡(nπ2)a_n = \int_{-0.5}^{0.5} \cos(n\pi t) dt = \frac{2}{n\pi} \sin\left( \frac{n\pi}{2} \right)an​=∫−0.50.5​cos(nπt)dt=nπ2​sin(2nπ​)（偶数项为0，奇数项交替为±2nπ\pm\frac{2}{n\pi}±nπ2​）；
• bn=0b_n = 0bn​=0（方波是偶函数，与奇函数sin⁡(nω0t)\sin(n\omega_0 t)sin(nω0​t)乘积的积分是0）。

因此，方波的傅里叶级数为：f(t)=12+2π(cos⁡(πt)−13cos⁡(3πt)+15cos⁡(5πt)−⋯ )f(t) = \frac{1}{2} + \frac{2}{\pi} \left( \cos(\pi t) - \frac{1}{3} \cos(3\pi t) + \frac{1}{5} \cos(5\pi t) - \cdots \right)f(t)=21​+π2​(cos(πt)−31​cos(3πt)+51​cos(5πt)−⋯)

当我们取前1项（基波）、前3项（基波+3次谐波）、前5项……时，级数的和会越来越接近原方波——这就是傅里叶级数的“逼近”效果！

5.4 复数形式：更简洁的数学表达

三角形式虽然直观，但计算繁琐。利用欧拉公式ejθ=cos⁡θ+jsin⁡θe^{j\theta} = \cos\theta + j\sin\thetaejθ=cosθ+jsinθ（j=−1j = \sqrt{-1}j=−1​是虚数单位），可以将傅里叶级数转化为复数形式：$
f(t) = \sum_{n=-\infty}^\infty c_n e^{j n \omega_0 t}
$其中cnc_ncn​是复数系数，表示第nnn次谐波的“复振幅”（包含振幅和相位信息）。

复数系数的计算公式为：cn=1T∫−T/2T/2f(t)e−jnω0tdt(n=0,±1,±2,… )c_n = \frac{1}{T} \int_{-T/2}^{T/2} f(t) e^{-j n \omega_0 t} dt \quad (n = 0, \pm1, \pm2, \dots)cn​=T1​∫−T/2T/2​f(t)e−jnω0​tdt(n=0,±1,±2,…)

复数形式的优势：

• 统一了余弦和正弦项，用单一复指数函数表示；
• 为傅里叶变换的推导提供了自然过渡；
• 频域表示更对称（包含正、负频率，但负频率是数学共轭，物理意义与正频率一致）。

6. 傅里叶变换：从周期到非周期的“飞跃”

傅里叶级数处理的是周期信号，但现实中的信号大多是非周期的（如一次闪电、一段语音）。如何扩展到非周期信号？

6.1 核心思路：非周期=“周期无穷大”

非周期信号可以看作“周期无穷大的周期信号”——当周期T→∞T \to \inftyT→∞时，信号不再重复，变成非周期。

6.2 从离散到连续：傅里叶变换的推导

假设我们有一个非周期函数f(t)f(t)f(t)，构造周期函数fT(t)f_T(t)fT​(t)（∣t∣<T/2|t| < T/2∣t∣<T/2时fT(t)=f(t)f_T(t) = f(t)fT​(t)=f(t)，否则重复）。当T→∞T \to \inftyT→∞时：

1. 频率离散→连续：角频率间隔Δω=ω0=2πT→0\Delta\omega = \omega_0 = \frac{2\pi}{T} \to 0Δω=ω0​=T2π​→0，离散频率nω0n\omega_0nω0​变成连续频率ω\omegaω；
2. 系数缩放：复数系数cn=1TF(nω0)c_n = \frac{1}{T} F(n\omega_0)cn​=T1​F(nω0​)，其中F(nω0)=∫−∞∞f(t)e−jnω0tdtF(n\omega_0) = \int_{-\infty}^\infty f(t) e^{-j n \omega_0 t} dtF(nω0​)=∫−∞∞​f(t)e−jnω0​tdt（当T→∞T \to \inftyT→∞时，F(nω0)→F(ω)F(n\omega_0) \to F(\omega)F(nω0​)→F(ω)，即傅里叶变换）。

6.3 傅里叶变换对：分解与重构

通过极限推导，最终得到傅里叶变换（正变换，时域→频域）和逆傅里叶变换（逆变换，频域→时域）的公式：

正变换（时域→频域）

F(ω)=F{f(t)}=∫−∞∞f(t)e−jωtdtF(\omega) = \mathcal{F}\{f(t)\} = \int_{-\infty}^\infty f(t) e^{-j\omega t} dtF(ω)=F{f(t)}=∫−∞∞​f(t)e−jωtdt

• F(ω)F(\omega)F(ω)称为f(t)f(t)f(t)的傅里叶变换（频域表示）；
• e−jωte^{-j\omega t}e−jωt是“基函数”（复指数波）；
• 积分表示将f(t)f(t)f(t)投影到所有频率ω\omegaω的基函数上，得到每个频率的“强度”。

逆变换（频域→时域）

f(t)=F−1{F(ω)}=12π∫−∞∞F(ω)ejωtdωf(t) = \mathcal{F}^{-1}\{F(\omega)\} = \frac{1}{2\pi} \int_{-\infty}^\infty F(\omega) e^{j\omega t} d\omegaf(t)=F−1{F(ω)}=2π1​∫−∞∞​F(ω)ejωtdω

• 逆变换是正变换的“逆操作”：将频域的所有复指数波叠加，重构出时域信号；
• 12π\frac{1}{2\pi}2π1​是归一化系数（补偿正变换的积分缩放）。

6.4 频域的物理意义：振幅谱与相位谱

傅里叶变换F(ω)F(\omega)F(ω)是复函数，可以表示为：F(ω)=∣F(ω)∣ejϕ(ω)F(\omega) = |F(\omega)| e^{j\phi(\omega)}F(ω)=∣F(ω)∣ejϕ(ω)其中：

• ∣F(ω)∣|F(\omega)|∣F(ω)∣：振幅谱（幅频特性），表示频率ω\omegaω的成分“有多强”；
• ϕ(ω)=arg⁡(F(ω))\phi(\omega) = \arg(F(\omega))ϕ(ω)=arg(F(ω))：相位谱（相频特性），表示该成分“在什么时候出现”。

6.5 关键性质：傅里叶的“实用技巧”

傅里叶变换有许多重要性质，以下是最常用的三个：

（1）线性性质

若F{f1(t)}=F1(ω)\mathcal{F}\{f_1(t)\} = F_1(\omega)F{f1​(t)}=F1​(ω)，F{f2(t)}=F2(ω)\mathcal{F}\{f_2(t)\} = F_2(\omega)F{f2​(t)}=F2​(ω)，则：F{af1(t)+bf2(t)}=aF1(ω)+bF2(ω)\mathcal{F}\{a f_1(t) + b f_2(t)\} = a F_1(\omega) + b F_2(\omega)F{af1​(t)+bf2​(t)}=aF1​(ω)+bF2​(ω)——复杂信号的傅里叶变换等于其成分的线性组合。

（2）时移性质

若F{f(t)}=F(ω)\mathcal{F}\{f(t)\} = F(\omega)F{f(t)}=F(ω)，则：F{f(t−t0)}=e−jωt0F(ω)\mathcal{F}\{f(t - t_0)\} = e^{-j\omega t_0} F(\omega)F{f(t−t0​)}=e−jωt0​F(ω)——时域信号延迟t0t_0t0​，振幅谱不变，相位谱增加−ωt0-\omega t_0−ωt0​（相位偏移与频率成正比）。

（3）频移性质

若F{f(t)}=F(ω)\mathcal{F}\{f(t)\} = F(\omega)F{f(t)}=F(ω)，则：F{f(t)ejω0t}=F(ω−ω0)\mathcal{F}\{f(t) e^{j\omega_0 t}\} = F(\omega - \omega_0)F{f(t)ejω0​t}=F(ω−ω0​)——时域信号乘以复指数ejω0te^{j\omega_0 t}ejω0​t，频域谱整体右移ω0\omega_0ω0​（通信中“调制”的理论基础）。

7. 适用边界与局限性

傅里叶变换是强大的工具，但并非“万能”——它有严格的适用条件和固有局限性。

7.1 适用条件：狄利克雷条件

一个函数f(t)f(t)f(t)能进行傅里叶变换，需满足狄利克雷条件：

1. 绝对可积：∫−∞∞∣f(t)∣dt<∞\int_{-\infty}^\infty |f(t)| dt < \infty∫−∞∞​∣f(t)∣dt<∞（充分非必要，如冲激函数δ(t)\delta(t)δ(t)不绝对可积，但傅里叶变换存在）；
2. 有限正则性：任意有限区间内，f(t)f(t)f(t)只有有限个极值点和有限个第一类不连续点（跳跃不连续，左、右极限存在）。

7.2 固有局限性：时间与频率的“权衡”

傅里叶变换是全局分析——它给出信号在整个时间域内的频率成分，但无法定位“某个频率在什么时候出现”。例如：

• 一个突然出现的高频脉冲，傅里叶变换会显示有高频成分，但无法确定脉冲的时间。

这源于海森堡不确定性原理：Δt⋅Δω≥12\Delta t \cdot \Delta\omega \geq \frac{1}{2}Δt⋅Δω≥21​

• Δt\Delta tΔt：时间分辨率（精确确定事件时间的能力）；
• Δω\Delta\omegaΔω：频率分辨率（精确区分频率成分的能力）；
• 结论：时间分辨率越高，频率分辨率越低；反之亦然。

7.3 不适用的场景

1. 非平稳信号：频率随时间变化的信号（如鸟鸣、音乐滑音），傅里叶变换无法捕捉频率的时变特性；
2. 强突变信号：含有尖锐突变的信号（如冲激的导数），傅里叶变换收敛慢；
3. 离散信号：连续傅里叶变换处理的是连续时间信号，离散信号需用离散傅里叶变换（DFT）或快速傅里叶变换（FFT）。

8. 傅里叶变换的“用武之地”

傅里叶变换是信号处理、图像处理、通信等领域的“基石”，以下是典型应用：

8.1 信号滤波：去掉噪声

例如，去掉语音中的高频噪声：

1. 时域信号f(t)=语音+高频噪声f(t) = \text{语音} + \text{高频噪声}f(t)=语音+高频噪声；
2. 正变换得到F(ω)=语音谱+噪声谱F(\omega) = \text{语音谱} + \text{噪声谱}F(ω)=语音谱+噪声谱；
3. 频域滤波：去掉高频成分（噪声）；
4. 逆变换得到干净的语音ffiltered(t)f_{\text{filtered}}(t)ffiltered​(t)。

8.2 图像处理：边缘检测

图像是二维信号（空间域的亮度分布），其傅里叶变换是二维傅里叶变换：F(u,v)=∫−∞∞∫−∞∞f(x,y)e−j(ux+vy)dxdyF(u,v) = \int_{-\infty}^\infty \int_{-\infty}^\infty f(x,y) e^{-j(ux + vy)} dx dyF(u,v)=∫−∞∞​∫−∞∞​f(x,y)e−j(ux+vy)dxdy

• 高频成分对应图像的边缘（亮度突变）；
• 提取高频成分，逆变换后得到边缘图像。

8.3 通信：调制与解调

低频信号（如语音）无法远距离传输，需调制到高频载波上：

1. 时域调制：s(t)=f(t)cos⁡(ωct)s(t) = f(t) \cos(\omega_c t)s(t)=f(t)cos(ωc​t)（ωc\omega_cωc​是载波频率）；
2. 频域效果：S(ω)=12F(ω−ωc)+12F(ω+ωc)S(\omega) = \frac{1}{2} F(\omega - \omega_c) + \frac{1}{2} F(\omega + \omega_c)S(ω)=21​F(ω−ωc​)+21​F(ω+ωc​)（低频谱复制到载波两侧）；
3. 接收端解调：去掉载波，恢复原信号（AM广播的原理）。

9. 总结：傅里叶变换的“本质”

傅里叶变换的核心是**“分解-重构”**的思维：

1. 分解：将时域的复杂信号拆成频域的正弦波（或复指数波）；
2. 处理：在频域对信号进行操作（如滤波、调制）；
3. 重构：逆变换回时域，得到处理后的信号。

对小白来说，记住这句话就够了：

傅里叶变换是“时域与频域之间的翻译器”——它把“时间的故事”翻译成“频率的乐谱”，让我们能从另一个角度理解世界。

最后，傅里叶变换不是终点，而是起点——它开启了“时频分析”的大门，后续的小波变换、短时傅里叶变换等，都是在傅里叶变换的基础上发展起来的。但无论工具多复杂，“分解与重构”的思想始终是核心。

附录：关键公式汇总

1. 傅里叶级数（三角形式）：f(t)=a02+∑n=1∞[ancos⁡(nω0t)+bnsin⁡(nω0t)]f(t) = \frac{a_0}{2} + \sum_{n=1}^\infty [a_n \cos(n\omega_0 t) + b_n \sin(n\omega_0 t)]f(t)=2a0​​+∑n=1∞​[an​cos(nω0​t)+bn​sin(nω0​t)]；
2. 傅里叶级数（复数形式）：f(t)=∑n=−∞∞cnejnω0tf(t) = \sum_{n=-\infty}^\infty c_n e^{j n \omega_0 t}f(t)=∑n=−∞∞​cn​ejnω0​t；
3. 傅里叶变换对：F(ω)=∫−∞∞f(t)e−jωtdt,f(t)=12π∫−∞∞F(ω)ejωtdωF(\omega) = \int_{-\infty}^\infty f(t) e^{-j\omega t} dt, \quad f(t) = \frac{1}{2\pi} \int_{-\infty}^\infty F(\omega) e^{j\omega t} d\omegaF(ω)=∫−∞∞​f(t)e−jωtdt,f(t)=2π1​∫−∞∞​F(ω)ejωtdω

希望这篇教程能帮你打开傅里叶变换的大门——接下来，就可以用它去解读更多“复杂世界的简单密码”了！

案例介绍

模拟一个由低频正弦波和高频噪声叠加的非周期信号，使用傅里叶变换分解信号得到振幅谱，通过频域滤波去除高频噪声后重构信号，对比滤波前后的时域波形差异。

# 导入必要的库importnumpyasnpimportmatplotlib.pyplotasplt# 定义生成模拟信号的函数defgenerate_signal(time:np.ndarray,f_low:float=5.0,# 低频信号频率，单位Hznoise_amp:float=0.3# 噪声振幅)->np.ndarray:""" 生成由低频正弦波和随机高频噪声叠加的模拟信号 :param time: 时间序列，一维数组 :param f_low: 低频正弦波频率，默认5Hz :param noise_amp: 高斯噪声振幅，默认0.3 :return: 叠加噪声后的模拟信号，一维数组 """# 生成低频正弦波low_freq_wave=np.sin(2*np.pi*f_low*time)# 生成随机高频噪声noise=noise_amp*np.random.randn(len(time))# 叠加得到模拟信号signal=low_freq_wave+noisereturnsignal# 定义傅里叶变换及滤波函数deffft_filter(signal:np.ndarray,sampling_rate:float,cutoff_freq:float=10.0# 滤波截止频率)->np.ndarray:""" 对信号进行傅里叶变换，频域滤波后重构信号 :param signal: 输入时域信号，一维数组 :param sampling_rate: 采样频率，单位Hz :param cutoff_freq: 低通滤波截止频率，默认10Hz :return: 滤波后的时域信号，一维数组 """# 对信号进行快速傅里叶变换fft_result=np.fft.fft(signal)# 获取频率轴freq_axis=np.fft.fftfreq(len(signal),d=1/sampling_rate)# 生成低通滤波掩码：保留截止频率以下的成分low_pass_mask=np.abs(freq_axis)<=cutoff_freq# 应用滤波掩码filtered_fft=fft_result*low_pass_mask# 逆傅里叶变换重构时域信号filtered_signal=np.real(np.fft.ifft(filtered_fft))returnfiltered_signal# 定义可视化函数defplot_results(time:np.ndarray,original:np.ndarray,filtered:np.ndarray,sampling_rate:float)->None:""" 绘制原始信号、滤波后信号及时域波形图 :param time: 时间序列，一维数组 :param original: 原始时域信号，一维数组 :param filtered: 滤波后时域信号，一维数组 :param sampling_rate: 采样频率，单位Hz :return: None """# 创建画布fig,axes=plt.subplots(2,1,figsize=(12,8))# 设置字体大小plt.rcParams['font.size']=12# 绘制原始信号axes[0].plot(time,original,color='blue',label='Original Signal')axes[0].set_title('Original Time-Domain Signal (with Noise)')axes[0].set_xlabel('Time (s)')axes[0].set_ylabel('Amplitude')axes[0].grid(True)axes[0].legend()# 绘制滤波后信号axes[1].plot(time,filtered,color='red',label='Filtered Signal')axes[1].set_title('Filtered Time-Domain Signal (without High-Frequency Noise)')axes[1].set_xlabel('Time (s)')axes[1].set_ylabel('Amplitude')axes[1].grid(True)axes[1].legend()# 调整子图间距plt.tight_layout()# 显示图像plt.show()# 主程序defmain():""" 主函数，执行信号生成、傅里叶滤波及结果可视化流程 :return: None """# 设置采样参数sampling_rate=100.0# 采样频率，单位Hzduration=2.0# 信号持续时间，单位s# 生成时间序列time=np.linspace(0.0,duration,int(sampling_rate*duration),endpoint=False)# 生成模拟信号original_signal=generate_signal(time)# 傅里叶变换频域滤波filtered_signal=fft_filter(original_signal,sampling_rate)# 可视化结果plot_results(time,original_signal,filtered_signal,sampling_rate)# 程序入口if__name__=="__main__":main()

代码功能总述

该代码实现了基于快速傅里叶变换(FFT)的频域低通滤波流程：

1. 生成“低频正弦波+随机高频高斯噪声”的混合模拟信号；
2. 对混合信号进行傅里叶变换得到频域谱；
3. 设计低通滤波器保留低频有用信号、滤除高频噪声；
4. 通过逆傅里叶变换重构时域信号，并可视化滤波前后的差异。

逐模块深度解析

1. 库导入模块

importnumpyasnpimportmatplotlib.pyplotasplt

• numpy (np)：提供向量/矩阵运算、傅里叶变换等数值计算核心功能，是科学计算的基础。
• matplotlib.pyplot (plt)：用于绘制时域信号波形图，可视化滤波效果。

2. 模拟信号生成函数generate_signal

defgenerate_signal(time:np.ndarray,f_low:float=5.0,# 低频信号频率，单位Hznoise_amp:float=0.3# 噪声振幅)->np.ndarray:# 生成低频正弦波low_freq_wave=np.sin(2*np.pi*f_low*time)# 生成随机高频噪声noise=noise_amp*np.random.randn(len(time))# 叠加得到模拟信号signal=low_freq_wave+noisereturnsignal

数学原理与代码逻辑：

• time：采样时间序列，由主函数通过np.linspace生成（等间隔采样点的时间戳）。
• 低频正弦波：公式为y(t)=sin⁡(2πflowt)y(t) = \sin(2\pi f_{\text{low}} t)y(t)=sin(2πflow​t)，其中2πflow2\pi f_{\text{low}}2πflow​是角频率，确保信号在单位时间内完成flowf_{\text{low}}flow​个周期。
• 高频噪声：用np.random.randn生成零均值高斯白噪声（高频特性源于其频谱平坦，包含所有频率成分），noise_amp控制噪声强度。
• 信号叠加：将有用低频波与噪声线性相加，模拟现实中“信号被噪声污染”的场景。

3. 傅里叶滤波核心函数fft_filter

deffft_filter(signal:np.ndarray,sampling_rate:float,cutoff_freq:float=10.0# 滤波截止频率)->np.ndarray:# 快速傅里叶变换：时域→频域fft_result=np.fft.fft(signal)# 生成频率轴freq_axis=np.fft.fftfreq(len(signal),d=1/sampling_rate)# 低通滤波掩码：保留截止频率以下的成分low_pass_mask=np.abs(freq_axis)<=cutoff_freq# 应用掩码滤波filtered_fft=fft_result*low_pass_mask# 逆傅里叶变换：频域→时域filtered_signal=np.real(np.fft.ifft(filtered_fft))returnfiltered_signal

这是代码的核心模块，涉及傅里叶变换的底层逻辑，需结合数学原理理解：

3.1 快速傅里叶变换(FFT)

• 函数：np.fft.fft(signal)
• 数学本质：将时域离散信号x(n)x(n)x(n)转换为频域离散谱X(k)X(k)X(k)，公式为：X(k)=∑n=0N−1x(n)e−j2πNknX(k) = \sum_{n=0}^{N-1} x(n) e^{-j\frac{2\pi}{N}kn}X(k)=∑n=0N−1​x(n)e−jN2π​kn
  其中NNN是信号长度，jjj是虚数单位。fft_result是复数数组，包含各频率成分的振幅和相位信息。

3.2 频率轴生成

• 函数：np.fft.fftfreq(len(signal), d=1/sampling_rate)
• 参数解释：
  • d：采样间隔，即相邻两个采样点的时间差Δt=1sampling_rate\Delta t = \frac{1}{\text{sampling\_rate}}Δt=sampling_rate1​；
  • 返回值：频率轴数组，范围为[−fs2,fs2)[-\frac{f_s}{2}, \frac{f_s}{2})[−2fs​​,2fs​​)（fsf_sfs​为采样频率），符合奈奎斯特采样定理（能处理的最高频率为fs2\frac{f_s}{2}2fs​​）。
• 例如：采样率fs=100Hzf_s=100\text{Hz}fs​=100Hz，则频率轴范围为[−50Hz,50Hz)[-50\text{Hz}, 50\text{Hz})[−50Hz,50Hz)。

3.3 低通滤波掩码

• 逻辑：low_pass_mask = np.abs(freq_axis) <= cutoff_freq
• 本质：生成一个布尔数组，仅保留绝对值≤截止频率的频率成分（即低频有用信号），其他高频噪声成分置零。
• 本例中：cutoff_freq=10\text{Hz}，因此保留[−10Hz,10Hz][-10\text{Hz}, 10\text{Hz}][−10Hz,10Hz]的频率，滤除高于10Hz的噪声（符合“低频5Hz信号+高频噪声”的场景设计）。

3.4 频域滤波与逆傅里叶变换

• 滤波：filtered_fft = fft_result * low_pass_mask（逐元素乘积，高频成分被掩码置零）；
• 逆FFT：np.fft.ifft(filtered_fft)将滤除后的频域信号转换回时域；
• 取实部：np.real(...)是因为浮点运算误差可能产生微小虚部，理论上时域信号应为实数，故需舍去虚部。

4. 可视化函数plot_results

defplot_results(time:np.ndarray,original:np.ndarray,filtered:np.ndarray,sampling_rate:float)->None:fig,axes=plt.subplots(2,1,figsize=(12,8))# 创建2行1列的画布plt.rcParams['font.size']=12# 设置字体大小# 绘制原始含噪信号axes[0].plot(time,original,color='blue',label='Original Signal')axes[0].set_title('Original Time-Domain Signal (with Noise)')axes[0].set_xlabel('Time (s)')axes[0].set_ylabel('Amplitude')axes[0].grid(True)axes[0].legend()# 绘制滤波后信号axes[1].plot(time,filtered,color='red',label='Filtered Signal')axes[1].set_title('Filtered Time-Domain Signal (without High-Frequency Noise)')axes[1].set_xlabel('Time (s)')axes[1].set_ylabel('Amplitude')axes[1].grid(True)axes[1].legend()plt.tight_layout()# 自动调整子图间距plt.show()# 显示图像

• 功能：将原始含噪信号与滤波后干净信号绘制在同一画布的上下两个子图中，直观对比滤波效果。
• 关键参数：figsize控制画布大小，grid=True添加网格线便于观察波形，label和legend用于标识曲线含义。

5. 主函数与程序入口

defmain():# 采样参数设置sampling_rate=100.0# 采样频率：每秒采集100个点duration=2.0# 信号持续时间：2秒# 生成时间序列：从0到2秒，共200个采样点（endpoint=False表示不含终点2.0）time=np.linspace(0.0,duration,int(sampling_rate*duration),endpoint=False)# 生成模拟含噪信号original_signal=generate_signal(time)# 傅里叶滤波filtered_signal=fft_filter(original_signal,sampling_rate)# 可视化结果plot_results(time,original_signal,filtered_signal,sampling_rate)if__name__=="__main__":main()

流程控制逻辑：

1. 采样参数：
   • 采样率fs=100Hzf_s=100\text{Hz}fs​=100Hz（满足奈奎斯特条件：有用信号频率5Hz，需fs>2×5=10Hzf_s > 2 \times 5 = 10\text{Hz}fs​>2×5=10Hz）；
   • 持续时间2秒，故采样点总数N=fs×duration=200N = f_s \times \text{duration} = 200N=fs​×duration=200。
2. 时间序列生成：np.linspace生成等间隔时间戳，endpoint=False确保时间点从0开始，步长为1/fs=0.01秒1/f_s=0.01\text{秒}1/fs​=0.01秒。
3. 函数调用链：生成信号 → 滤波 → 可视化，构成完整的“信号处理流水线”。

关键数学建模知识点总结

运行效果与结论

代码运行后将输出两张时域波形图：

• 上方：原始含噪信号（波形受高频噪声干扰，呈现“毛刺”）；
• 下方：滤波后信号（高频噪声被滤除，恢复为平滑的5Hz正弦波）。

该流程展示了傅里叶变换在信号去噪中的经典应用，是数学建模中“信号处理类问题”的基础方法之一。