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

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

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(ω0t) 的周期 T=2πω0T = \frac{2\pi}{\omega_0}T=ω02π
  • 其中 ω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)abg(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(ω0t),sin(ω0t),cos(2ω0t),sin(2ω0t),}这个函数系在 [−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 = 0T/2T/2cos(mω0t)sin(nω0t)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 = 0T/2T/2cos(mω0t)cos(nω0t)dt=0meqnm 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[ancos(nω0t)+bnsin(nω0t)]其中:

  • a02\frac{a_0}{2}2a0直流分量(常数项,对应频率0的成分);
  • ancos⁡(nω0t)a_n \cos(n\omega_0 t)ancos(nω0t)余弦谐波bnsin⁡(nω0t)b_n \sin(n\omega_0 t)bnsin(nω0t)正弦谐波
  • 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/2f(t)dt=T/2T/22a0dt+n=1[anT/2T/2cos(nω0t)dt+bnT/2T/2sin(nω0t)dt]由于余弦和正弦在整周期内的积分是0,右边求和项全为0,因此:a0=2T∫−T/2T/2f(t)dta_0 = \frac{2}{T} \int_{-T/2}^{T/2} f(t) dta0=T2T/2T/2f(t)dt

(2)计算余弦系数 ana_nan

将傅里叶级数两边乘以 cos⁡(nω0t)\cos(n\omega_0 t)cos(nω0t) 并积分:∫−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/2f(t)cos(nω0t)dt=an2T(其他项因正交性消失)因此: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=T2T/2T/2f(t)cos(nω0t)dt(n1)

(3)计算正弦系数 bnb_nbn

类似地,乘以 sin⁡(nω0t)\sin(n\omega_0 t)sin(nω0t) 并积分: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=T2T/2T/2f(t)sin(nω0t)dt(n1)

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)={10t<0.50.5<t<1

计算系数:

  • a0=∫−11f(t)dt=1a_0 = \int_{-1}^{1} f(t) dt = 1a0=11f(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.5cos(t)dt=2sin(2)(偶数项为0,奇数项交替为 ±2nπ\pm\frac{2}{n\pi}±2);
  • bn=0b_n = 0bn=0(方波是偶函数,与奇函数 sin⁡(nω0t)\sin(n\omega_0 t)sin(nω0t) 乘积的积分是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)31cos(3πt)+51cos(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=T1T/2T/2f(t)ejnω0tdt(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/2t<T/2fT(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=T1F(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)ejnω0tdt(当 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)etdt

  • F(ω)F(\omega)F(ω) 称为 f(t)f(t)f(t)傅里叶变换(频域表示);
  • e−jωte^{-j\omega t}et 是“基函数”(复指数波);
  • 积分表示将 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)=F1{F(ω)}=2π1F(ω)etdω

  • 逆变换是正变换的“逆操作”:将频域的所有复指数波叠加,重构出时域信号;
  • 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(tt0)}=et0F(ω)——时域信号延迟 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ω0t}=F(ωω0)——时域信号乘以复指数 ejω0te^{j\omega_0 t}ejω0t,频域谱整体右移 ω0\omega_0ω0(通信中“调制”的理论基础)。

7. 适用边界与局限性

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

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

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

  1. 绝对可积∫−∞∞∣f(t)∣dt<∞\int_{-\infty}^\infty |f(t)| dt < \inftyf(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)ej(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(ωct)ω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(ω)=21F(ωωc)+21F(ω+ω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[ancos(nω0t)+bnsin(nω0t)]
  2. 傅里叶级数(复数形式):f(t)=∑n=−∞∞cnejnω0tf(t) = \sum_{n=-\infty}^\infty c_n e^{j n \omega_0 t}f(t)=n=cnejnω0t
  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)etdt,f(t)=2π1F(ω)etdω

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

案例介绍

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

# 导入必要的库
import numpy as np
import matplotlib.pyplot as plt

# 定义生成模拟信号的函数
def generate_signal(
    time: np.ndarray,
    f_low: float = 5.0,  # 低频信号频率,单位Hz
    noise_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 + noise
    return signal

# 定义傅里叶变换及滤波函数
def fft_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))
    return filtered_signal

# 定义可视化函数
def plot_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()

# 主程序
def main():
    """
    主函数,执行信号生成、傅里叶滤波及结果可视化流程
    :return: None
    """
    # 设置采样参数
    sampling_rate = 100.0  # 采样频率,单位Hz
    duration = 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. 库导入模块
import numpy as np
import matplotlib.pyplot as plt
  • numpy (np):提供向量/矩阵运算、傅里叶变换等数值计算核心功能,是科学计算的基础。
  • matplotlib.pyplot (plt):用于绘制时域信号波形图,可视化滤波效果。
2. 模拟信号生成函数 generate_signal
def generate_signal(
    time: np.ndarray,
    f_low: float = 5.0,  # 低频信号频率,单位Hz
    noise_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 + noise
    return signal

数学原理与代码逻辑

  • time:采样时间序列,由主函数通过np.linspace生成(等间隔采样点的时间戳)。
  • 低频正弦波:公式为 y(t)=sin⁡(2πflowt)y(t) = \sin(2\pi f_{\text{low}} t)y(t)=sin(2πflowt),其中 2πflow2\pi f_{\text{low}}2πflow 是角频率,确保信号在单位时间内完成 flowf_{\text{low}}flow 个周期。
  • 高频噪声:用np.random.randn生成零均值高斯白噪声(高频特性源于其频谱平坦,包含所有频率成分),noise_amp控制噪声强度。
  • 信号叠加:将有用低频波与噪声线性相加,模拟现实中“信号被噪声污染”的场景。
3. 傅里叶滤波核心函数 fft_filter
def fft_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))
    return filtered_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=0N1x(n)ejN2π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(逐元素乘积,高频成分被掩码置零);
  • 逆FFTnp.fft.ifft(filtered_fft) 将滤除后的频域信号转换回时域;
  • 取实部np.real(...) 是因为浮点运算误差可能产生微小虚部,理论上时域信号应为实数,故需舍去虚部。
4. 可视化函数 plot_results
def plot_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 添加网格线便于观察波形,labellegend 用于标识曲线含义。
5. 主函数与程序入口
def main():
    # 采样参数设置
    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. 函数调用链:生成信号 → 滤波 → 可视化,构成完整的“信号处理流水线”。

关键数学建模知识点总结

代码模块核心数学原理建模意义
信号生成正弦波公式 $\sin(2\pi ft)$、高斯白噪声特性模拟现实中“低信噪比”的非周期信号场景
FFT变换离散傅里叶变换(DFT)、奈奎斯特采样定理将时域“混叠信号”分解为频域“频谱成分”,为滤波提供基础
频域滤波时域卷积=频域乘积、低通滤波器设计直接在频域分离有用信号与噪声,比时域滤波更高效直观
逆FFT傅里叶逆变换(IDFT)重构滤波后的干净时域信号

运行效果与结论

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

  • 上方:原始含噪信号(波形受高频噪声干扰,呈现“毛刺”);
  • 下方:滤波后信号(高频噪声被滤除,恢复为平滑的5Hz正弦波)。

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

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

您好啊数模君

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值