精华内容
下载资源
问答
  • 实验傅里叶变换系统的频域分析,实验代码,实验报告
  • 信号与系统实验仿真系统的MATLAB实现-傅里叶变换曲线(振幅频谱).m 《信号与系统实验仿真系统的MATLAB实现 程序的名字说明了程序的功能:)
  • 傅里叶变换、DFT、FFT分析理解;

    一、前言

    CSDN手机上查看数学公式,分式容易重叠,可以查看博客园文章:[信号与系统]傅里叶变换、DFT、FFT分析与理解

    说来惭愧,本科专业是测控技术与仪器,四舍五入也算个电子信息人(虽然我不怎么承认,毕竟本科也没怎么听课,天天在实验室玩单片机,实在是给专业丢人)。
    前不久因为工作原因(实际上是因为接别人的锅,别人忙,所以就让我这个菜鸡接人家剩下的继续搞,两个周后出程序)。用户要求输出电压、电流周波数据,说是用于负荷分析。别的我不用关心,我只要能输出周波数据就行。但公司有流程,所有的程序出厂前都需要进行测试、出厂检等等步骤,考虑到测试组能力有限,所以给他们做了一个软件,模拟用户对数据进行处理(实际上就是简单的进行傅里叶变换,得到频谱然后显示,测试组看个幅值就可以了)。但最终软件做完了,测试组说只有138台,他们没时间测。好吧,写一篇博客记录一下整个学习过程。

    二、傅里叶变换

    1.傅里叶级数

    最早接触到的傅里叶变换,估计大家都一样,都是从大学高数书上的,泰勒展开、级数后面就是这个玩意。我但是学的书,这章是打※号的,说白了就是不怎么重要,而且讲高数的老师说你们学信息的以后专业课上会详细学到(这句话好熟悉),但是只记得怎么解题,总结来说就是下面的这个公式:

    公式1: f ( t ) = A 0 + ∑ n = 1 ∞ [ A n ⋅ s i n ( n ω t + ϕ n ) ] f(t) = A_0 + \sum _{n=1}^{\infin} [A_n \cdot sin(n\omega t + \phi_n)] f(t)=A0+n=1[Ansin(nωt+ϕn)]

    更多的时候,我们见到的是下面的式子:

    公式2: f ( t ) = a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n π t l ) + b n ⋅ s i n ( n π t l ) ] f(t) = \frac{a_0}{2} + \sum _{n=1}^\infin [a_n \cdot cos(\frac{n\pi t}{l}) + b_n \cdot sin(\frac{n\pi t}{l})] f(t)=2a0+n=1[ancos(lnπt)+bnsin(lnπt)]

    π t l = x \frac{\pi t}{l} = x lπt=x即可得到:

    公式3: f ( t ) = a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n x ) + b n ⋅ s i n ( n x ) ] f(t) = \frac{a_0}{2} + \sum _{n=1}^\infin [a_n \cdot cos(nx) + b_n \cdot sin(nx)] f(t)=2a0+n=1[ancos(nx)+bnsin(nx)]

    由公式1到公式2的转换,简单证明:

    利用公式: s i n ( a + b ) = s i n ( a ) c o s ( b ) + c o s ( a ) s i n ( b ) sin(a+b) = sin(a)cos(b) + cos(a)sin(b) sin(a+b)=sin(a)cos(b)+cos(a)sin(b)

    得到: A n ⋅ s i n ( n ω t + ϕ n ) = A n ⋅ s i n ( ϕ n ) ⋅ c o s ( n ω t ) + A n ⋅ c o s ( ϕ n ) ⋅ s i n ( n ω t ) A_n \cdot sin(n\omega t + \phi_n) = A_n \cdot sin(\phi _n)\cdot cos(n\omega t) + A_n \cdot cos(\phi _n) \cdot sin(n\omega t) Ansin(nωt+ϕn)=Ansin(ϕn)cos(nωt)+Ancos(ϕn)sin(nωt)

    a 0 2 = A 0 , a n = A n ⋅ s i n ( ϕ n ) , b n = A n ⋅ c o s ( ϕ n ) , ω = π l \frac{a_0}{2} = A_0,a_n = A_n \cdot sin(\phi _n),b_n = A_n \cdot cos(\phi _n),\omega = \frac{\pi}{l} 2a0=A0an=Ansin(ϕn)bn=Ancos(ϕn)ω=lπ

    即可得到: f ( t ) = a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n π t l ) + b n ⋅ s i n ( n π t l ) ] f(t) = \frac{a_0}{2} + \sum _{n=1}^\infin [a_n \cdot cos(\frac{n\pi t}{l}) + b_n \cdot sin(\frac{n\pi t}{l})] f(t)=2a0+n=1[ancos(lnπt)+bnsin(lnπt)]

    π t l = x \frac{\pi t}{l} = x lπt=x,得到: f ( t ) = a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n x ) + b n ⋅ s i n ( n x ) ] f(t) = \frac{a_0}{2} + \sum _{n=1}^\infin [a_n \cdot cos(nx) + b_n \cdot sin(nx)] f(t)=2a0+n=1[ancos(nx)+bnsin(nx)]

    简单的说,就是任何一个周期信号都可被一系列三角函数表示。(至于为什么,我也不太清楚,暂且当回拿来主义吧,参考知乎大神:傅里叶系列(一)傅里叶级数的推导

    2.傅里叶级数系数求解

    对于傅里叶级数,我们要求的就是这些系数: a n 、 b n a_n、b_n anbn

    2.1.求解方法

    解法也是很简单,利用三角函数的正交性,对公式3两边积分
    ∫ − ∞ + ∞ f ( x ) ⋅ s i n ( n x ) d x = ∫ − ∞ + ∞ { a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n x ) + b n ⋅ s i n ( n x ) ] } ⋅ s i n ( n x ) d x \int_{-\infin}^{+\infin}f(x)\cdot sin(nx)dx = \int_{-\infin}^{+\infin}\{\frac{a_0}{2} + \sum_{n=1}^\infin [a_n \cdot cos(nx) + b_n \cdot sin(nx)]\}\cdot sin(nx)dx +f(x)sin(nx)dx=+{2a0+n=1[ancos(nx)+bnsin(nx)]}sin(nx)dx

    ∫ − ∞ + ∞ f ( x ) ⋅ c o s ( n x ) d x = ∫ − ∞ + ∞ { a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n x ) + b n ⋅ s i n ( n x ) ] } ⋅ c o s ( n x ) d x \int_{-\infin}^{+\infin}f(x)\cdot cos(nx)dx = \int_{-\infin}^{+\infin}\{\frac{a_0}{2} + \sum_{n=1}^\infin [a_n \cdot cos(nx) + b_n \cdot sin(nx)]\}\cdot cos(nx)dx +f(x)cos(nx)dx=+{2a0+n=1[ancos(nx)+bnsin(nx)]}cos(nx)dx

    2.2.三角函数的正交性

    利用三角函数的正交性:
    所谓正交性,书本概念请自行搜索,我只说简单提一句:正交性就是相乘为0

    ∫ − T 2 T 2 s i n ( n x ) ⋅ f ( k x ) d x = 0 \int_{-\frac{T}{2}}^{\frac{T}{2}}sin(nx)\cdot f(kx) dx = 0 2T2Tsin(nx)f(kx)dx=0
    当 f ( k x ) = 0 , 1 , s i n ( 1 x ) , c o s ( 1 x ) , s i n ( 2 x ) , c o s ( 2 x ) , … , s i n ( n x ) , c o s ( n x ) , 且 n ≠ k 当f(kx) = 0,1,sin(1x),cos(1x),sin(2x),cos(2x),\dots,sin(nx),cos(nx),且n\neq k f(kx)=0,1,sin(1x),cos(1x),sin(2x),cos(2x),,sin(nx),cos(nx)n=k
    其 中 : 0 = s i x ( 0 ⋅ x ) , 1 = c o s ( 0 ⋅ x ) 其中:0=six(0\cdot x),1=cos(0\cdot x) 0=six(0x)1=cos(0x)

    ∫ − T 2 T 2 c o s ( n x ) ⋅ f ( k x ) d x = 0 \int_{-\frac{T}{2}}^{\frac{T}{2}}cos(nx)\cdot f(kx) dx = 0 2T2Tcos(nx)f(kx)dx=0
    当 f ( k x ) = 0 , 1 , s i n ( 1 x ) , c o s ( 1 x ) , s i n ( 2 x ) , c o s ( 2 x ) , … , s i n ( n x ) , c o s ( n x ) , 且 n ≠ k 当f(kx) = 0,1,sin(1x),cos(1x),sin(2x),cos(2x),\dots,sin(nx),cos(nx),且n\neq k f(kx)=0,1,sin(1x),cos(1x),sin(2x),cos(2x),,sin(nx),cos(nx)n=k
    其 中 : 0 = s i x ( 0 ⋅ x ) , 1 = c o s ( 0 ⋅ x ) 其中:0=six(0\cdot x),1=cos(0\cdot x) 0=six(0x)1=cos(0x)

    即,当 n = k n = k n=k时, ∫ − T 2 T 2 c o s ( n x ) ⋅ f ( n x ) d x ≠ 0 \int_{-\frac{T}{2}}^{\frac{T}{2}}cos(nx)\cdot f(nx) dx \neq 0 2T2Tcos(nx)f(nx)dx=0

    关于这个正交性很好证明,积分+倍角变换很容易就出来了,我这里就不去证明了。

    2.3.系数求解过程

    接下来我主要解释以下如何得到系数:

    首先,求 a 0 a_0 a0

    ∫ − ∞ + ∞ f ( x ) ⋅ c o s ( 0 ) d x = ∫ − ∞ + ∞ { a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n x ) + b n ⋅ s i n ( n x ) ] } ⋅ c o s ( 0 ) d x \int_{-\infin}^{+\infin}f(x)\cdot cos(0)dx = \int_{-\infin}^{+\infin}\{\frac{a_0}{2} + \sum_{n=1}^\infin [a_n \cdot cos(nx) + b_n \cdot sin(nx)]\}\cdot cos(0)dx +f(x)cos(0)dx=+{2a0+n=1[ancos(nx)+bnsin(nx)]}cos(0)dx

    由于 f ( x ) f(x) f(x)本身就是周期的,加上 c o s ( n x ) cos(nx) cos(nx)也是周期的,因此只需要求一个周期内的即可:

    ∫ − T 2 T 2 f ( x ) ⋅ c o s ( 0 ) d x = ∫ − T 2 T 2 { a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n x ) + b n ⋅ s i n ( n x ) ] } ⋅ c o s ( 0 ) d x \int_{-\frac{T}{2}}^{\frac{T}{2}}f(x)\cdot cos(0)dx = \int_{-\frac{T}{2}}^{\frac{T}{2}}\{\frac{a_0}{2} + \sum_{n=1}^\infin [a_n \cdot cos(nx) + b_n \cdot sin(nx)]\}\cdot cos(0)dx 2T2Tf(x)cos(0)dx=2T2T{2a0+n=1[ancos(nx)+bnsin(nx)]}cos(0)dx

    利用三角函数的正交性,以及 c o s ( 0 ) = 1 cos(0)=1 cos(0)=1,得:

    ∫ − T 2 T 2 f ( x ) ⋅ 1 d x = ∫ − T 2 T 2 a 0 2 ⋅ 1 d x \int_{-\frac{T}{2}}^{\frac{T}{2}}f(x)\cdot 1dx = \int_{-\frac{T}{2}}^{\frac{T}{2}}\frac{a_0}{2} \cdot 1dx 2T2Tf(x)1dx=2T2T2a01dx

    这里的 T = 2 π T=2\pi T=2π,然后因为 a 0 2 \frac{a_0}{2} 2a0不是被积量,所以替换、系数提出:

    ∫ − π π f ( x ) ⋅ 1 d x = a 0 2 ⋅ ∫ − π π 1 d x \int_{-\pi}^{\pi}f(x)\cdot 1dx = \frac{a_0}{2} \cdot \int_{-\pi}^{\pi} 1dx ππf(x)1dx=2a0ππ1dx

    然后就一目了然了:

    ∫ − π π 1 d x = 2 π \int_{-\pi}^{\pi} 1dx = 2\pi ππ1dx=2π

    所以:

    a 0 = 1 π ⋅ ∫ − π π f ( x ) d x a_0 = \frac{1}{\pi} \cdot \int_{-\pi}^{\pi}f(x)dx a0=π1ππf(x)dx

    然后是计算 a n a_n an b n b_n bn
    同理,只不过带入的是 n n n

    ∫ − T 2 T 2 f ( x ) ⋅ c o s ( n x ) d x = ∫ − T 2 T 2 { a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n x ) + b n ⋅ s i n ( n x ) ] } ⋅ c o s ( n x ) d x \int_{-\frac{T}{2}}^{\frac{T}{2}}f(x)\cdot cos(nx)dx = \int_{-\frac{T}{2}}^{\frac{T}{2}}\{\frac{a_0}{2} + \sum_{n=1}^\infin [a_n \cdot cos(nx) + b_n \cdot sin(nx)]\}\cdot cos(nx)dx 2T2Tf(x)cos(nx)dx=2T2T{2a0+n=1[ancos(nx)+bnsin(nx)]}cos(nx)dx

    同样利用的三角函数的正交性,得:

    ∫ − T 2 T 2 f ( x ) ⋅ c o s ( n x ) d x = ∫ − T 2 T 2 [ a 0 2 + a n ⋅ c o s ( n x ) ] ⋅ c o s ( n x ) d x \int_{-\frac{T}{2}}^{\frac{T}{2}}f(x)\cdot cos(nx)dx = \int_{-\frac{T}{2}}^{\frac{T}{2}}[\frac{a_0}{2} +a_n\cdot cos(nx)]\cdot cos(nx)dx 2T2Tf(x)cos(nx)dx=2T2T[2a0+ancos(nx)]cos(nx)dx

    实际上应该这么写:

    ∫ − T 2 T 2 f ( x ) ⋅ c o s ( n x ) d x = ∫ − T 2 T 2 [ a 0 2 ⋅ 1 + a n ⋅ c o s ( n x ) ] ⋅ c o s ( n x ) d x \int_{-\frac{T}{2}}^{\frac{T}{2}}f(x)\cdot cos(nx)dx = \int_{-\frac{T}{2}}^{\frac{T}{2}}[\frac{a_0}{2} \cdot 1+a_n\cdot cos(nx)]\cdot cos(nx)dx 2T2Tf(x)cos(nx)dx=2T2T[2a01+ancos(nx)]cos(nx)dx

    发现漏掉了一个正交为0的项:

    ∫ − T 2 T 2 a 0 2 ⋅ [ 1 ⋅ c o s ( n x ) ] d x = a 0 2 ⋅ ∫ − T 2 T 2 [ 1 ⋅ c o s ( n x ) ] d x = 0 \int_{-\frac{T}{2}}^{\frac{T}{2}}\frac{a_0}{2} \cdot [1\cdot cos(nx)]dx = \frac{a_0}{2} \cdot \int_{-\frac{T}{2}}^{\frac{T}{2}}[1\cdot cos(nx)]dx = 0 2T2T2a0[1cos(nx)]dx=2a02T2T[1cos(nx)]dx=0

    所以,就简化为:

    ∫ − T 2 T 2 f ( x ) ⋅ c o s ( n x ) d x = ∫ − T 2 T 2 [ a n ⋅ c o s ( n x ) ] ⋅ c o s ( n x ) d x \int_{-\frac{T}{2}}^{\frac{T}{2}}f(x)\cdot cos(nx)dx = \int_{-\frac{T}{2}}^{\frac{T}{2}}[a_n\cdot cos(nx)]\cdot cos(nx)dx 2T2Tf(x)cos(nx)dx=2T2T[ancos(nx)]cos(nx)dx

    然后,按照和求 a 0 a_0 a0一样的处理:提出 a n a_n an、再积分…

    ∫ − T 2 T 2 f ( x ) ⋅ c o s ( n x ) d x = a n ⋅ ∫ − T 2 T 2 [ c o s ( n x ) ⋅ c o s ( n x ) ] d x \int_{-\frac{T}{2}}^{\frac{T}{2}}f(x)\cdot cos(nx)dx = a_n\cdot \int_{-\frac{T}{2}}^{\frac{T}{2}}[cos(nx)\cdot cos(nx)]dx 2T2Tf(x)cos(nx)dx=an2T2T[cos(nx)cos(nx)]dx

    利用倍角公式: 1 + c o s ( 2 x ) = 2 c o s 2 ( x ) 1+cos(2x)=2cos^2(x) 1+cos(2x)=2cos2(x)

    ∫ − T 2 T 2 [ c o s ( n x ) ⋅ c o s ( n x ) ] d x = ∫ − T 2 T 2 { 1 2 ⋅ [ 1 + c o s ( 2 n x ) ] } d x \int_{-\frac{T}{2}}^{\frac{T}{2}}[cos(nx)\cdot cos(nx)]dx = \int_{-\frac{T}{2}}^{\frac{T}{2}}\{\frac{1}{2}\cdot[1+cos(2nx)]\}dx 2T2T[cos(nx)cos(nx)]dx=2T2T{21[1+cos(2nx)]}dx

    然后还是利用正交性: ∫ − T 2 T 2 [ 1 ⋅ c o s ( 2 n x ) ] d x = 0 \int_{-\frac{T}{2}}^{\frac{T}{2}}[ 1 \cdot cos(2nx)]dx = 0 2T2T[1cos(2nx)]dx=0

    啰嗦一句:
    ∫ − T 2 T 2 [ 1 ⋅ c o s ( 2 n x ) ] d x = 0 \int_{-\frac{T}{2}}^{\frac{T}{2}}[ 1 \cdot cos(2nx)]dx = 0 2T2T[1cos(2nx)]dx=0式中相当于把 c o s ( n x ) cos(nx) cos(nx)倍频了,即一个周期T内塞下了两倍的 c o s ( n x ) cos(nx) cos(nx)信号
    因为 1 ⋅ c o s ( n x ) 1\cdot cos(nx) 1cos(nx)在一个周期内积分 = 0
    所以 1 ⋅ c o s ( 2 n x ) 1\cdot cos(2nx) 1cos(2nx)在一个周期内积分也 = 0

    所以,就剩下了: ∫ − T 2 T 2 f ( x ) ⋅ c o s ( n x ) d x = ∫ − T 2 T 2 [ 1 2 ⋅ a n ] d x \int_{-\frac{T}{2}}^{\frac{T}{2}}f(x)\cdot cos(nx)dx = \int_{-\frac{T}{2}}^{\frac{T}{2}}[\frac{1}{2}\cdot a_n]dx 2T2Tf(x)cos(nx)dx=2T2T[21an]dx

    整理得到: a n = 2 T ⋅ ∫ − T 2 T 2 f ( x ) ⋅ c o s ( n x ) d x a_n =\frac{2}{T} \cdot \int_{-\frac{T}{2}}^{\frac{T}{2}}f(x)\cdot cos(nx)dx an=T22T2Tf(x)cos(nx)dx

    同样 T = 2 π T=2\pi T=2π,所以替换:

    a n = 2 2 π ⋅ ∫ − π π f ( x ) ⋅ c o s ( n x ) d x = 1 π ⋅ ∫ − π π f ( x ) ⋅ c o s ( n x ) d x a_n = \frac{2}{2\pi} \cdot \int_{-\pi}^{\pi}f(x)\cdot cos(nx)dx = \frac{1}{\pi} \cdot \int_{-\pi}^{\pi}f(x)\cdot cos(nx)dx an=2π2ππf(x)cos(nx)dx=π1ππf(x)cos(nx)dx

    再把之前算出来的 a 0 a_0 a0拿过来:

    a 0 = 1 π ⋅ ∫ − π π f ( x ) d x a_0 = \frac{1}{\pi} \cdot \int_{-\pi}^{\pi}f(x)dx a0=π1ππf(x)dx

    a n = 1 π ⋅ ∫ − π π f ( x ) ⋅ c o s ( n x ) d x a_n = \frac{1}{\pi} \cdot \int_{-\pi}^{\pi}f(x)\cdot cos(nx)dx an=π1ππf(x)cos(nx)dx

    其实可以统一用: a n = 1 π ⋅ ∫ − π π f ( x ) ⋅ c o s ( n x ) d x a_n = \frac{1}{\pi} \cdot \int_{-\pi}^{\pi}f(x)\cdot cos(nx)dx an=π1ππf(x)cos(nx)dx表示

    同理: b n = 1 π ⋅ ∫ − π π f ( x ) ⋅ s i n ( n x ) d x b_n = \frac{1}{\pi} \cdot \int_{-\pi}^{\pi}f(x)\cdot sin(nx)dx bn=π1ππf(x)sin(nx)dx

    2.4.关于傅里叶级数的个人感悟

    但是,很多教材上往往都是令这个等于那个,令那个等于这个,弄得云里雾里,比如前面我列举了:
    傅里叶老先生费劲千辛万苦想出来: f ( x ) = A 0 + ∑ n = 1 ∞ [ A n ⋅ s i n ( n x + ϕ n ) ] f(x) = A_0 + \sum_{n=1}^{\infin} [A_n \cdot sin(nx + \phi_n)] f(x)=A0+n=1[Ansin(nx+ϕn)]
    后来的数学大师为了方便理解,转化出: f ( x ) = a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n x ) + b n ⋅ s i n ( n x ) ] f(x) = \frac{a_0}{2} + \sum_{n=1}^\infin [a_n \cdot cos(nx) + b_n \cdot sin(nx)] f(x)=2a0+n=1[ancos(nx)+bnsin(nx)]
    这两个式子我在开篇就写过,就是通过三角变换得到:把表示相位的分量 ϕ n \phi _n ϕn去掉了,转而用两个基本的正交的三角函数—— c o s ( x ) cos(x) cos(x) s i n ( x ) sin(x) sin(x)——表示。

    但是,为什么用 a 0 2 \frac{a_0}{2} 2a0而不是 a 0 a_0 a0
    本身 a 0 a_0 a0就是一个待求系数, a 0 2 \frac{a_0}{2} 2a0或者 a 0 a_0 a0,你用哪一个都一样啊。

    反正我本科学高数的时候,后来学信号与处理的时候,甚至再到后面考研刷数学,我都没有去想过这个。直到工作后重新学这个傅里叶变换,才真正静下心来去思考这个问题。

    首先列举两个傅里叶级数的基本公式:

    式子1: f ( x ) = a 0 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n x ) + b n ⋅ s i n ( n x ) ] f(x) = a_0 + \sum _{n=1}^\infin [a_n \cdot cos(nx) + b_n \cdot sin(nx)] f(x)=a0+n=1[ancos(nx)+bnsin(nx)]

    式子2: f ( x ) = a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n x ) + b n ⋅ s i n ( n x ) ] f(x) = \frac{a_0}{2} + \sum _{n=1}^\infin [a_n \cdot cos(nx) + b_n \cdot sin(nx)] f(x)=2a0+n=1[ancos(nx)+bnsin(nx)]

    可以看出来,两个式子的区别就是式子1的第一项是 a 0 a_0 a0,式子2的第一项是 a 0 2 \frac{a_0}{2} 2a0

    然后按照前面我们求系数的方式来求解这个 a 0 a_0 a0注意,两个式子都是求 a 0 a_0 a0

    式子1: a 0 ⋅ ∫ − T 2 T 2 d x = ∫ − T 2 T 2 f ( x ) d x a_0\cdot \int _{-\frac{T}{2}}^{\frac{T}{2}}dx = \int _{-\frac{T}{2}}^{\frac{T}{2}}f(x)dx a02T2Tdx=2T2Tf(x)dx

    式子2: a 0 2 ⋅ ∫ − T 2 T 2 d x = ∫ − T 2 T 2 f ( x ) d x \frac{a_0}{2}\cdot \int _{-\frac{T}{2}}^{\frac{T}{2}}dx = \int _{-\frac{T}{2}}^{\frac{T}{2}}f(x)dx 2a02T2Tdx=2T2Tf(x)dx

    这里就不去对这个周期 T T T进行处理了,就用 T T T这个符号来表示:

    式子1: a 0 = 1 T ⋅ ∫ − T 2 T 2 f ( x ) d x a_0 = \frac{1}{T}\cdot \int _{-\frac{T}{2}}^{\frac{T}{2}}f(x)dx a0=T12T2Tf(x)dx

    式子2: a 0 = 2 T ⋅ ∫ − T 2 T 2 f ( x ) d x a_0 = \frac{2}{T}\cdot \int _{-\frac{T}{2}}^{\frac{T}{2}}f(x)dx a0=T22T2Tf(x)dx

    然后对比之前求的 a n a_n an a n = 2 T ⋅ ∫ − T 2 T 2 f ( x ) ⋅ c o s ( n x ) d x a_n = \frac{2}{T} \cdot \int _{-\frac{T}{2}}^{\frac{T}{2}}f(x)\cdot cos(nx)dx an=T22T2Tf(x)cos(nx)dx

    会发现:式子2中那种用 a 0 2 \frac{a_0}{2} 2a0表示的系数,求的 a 0 a_0 a0可以直接套用 a n a_n an的求解公式求出,得到的这个 a 0 a_0 a0并不是傅里叶公式中的系数,你还需要除以2,因为给你的傅里叶表达式是 f ( x ) = a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n x ) + b n ⋅ s i n ( n x ) ] f(x) = \frac{a_0}{2} + \sum _{n=1}^\infin [a_n \cdot cos(nx) + b_n \cdot sin(nx)] f(x)=2a0+n=1[ancos(nx)+bnsin(nx)]

    总之用 f ( x ) = a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n x ) + b n ⋅ s i n ( n x ) ] f(x) = \frac{a_0}{2} + \sum _{n=1}^\infin [a_n \cdot cos(nx) + b_n \cdot sin(nx)] f(x)=2a0+n=1[ancos(nx)+bnsin(nx)],你只需要记住两个式子(如果是需要应付考试速成的,请直接用下面的三个公式):

    ①: f ( x ) = a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n x ) + b n ⋅ s i n ( n x ) ] f(x) = \frac{a_0}{2} + \sum _{n=1}^\infin [a_n \cdot cos(nx) + b_n \cdot sin(nx)] f(x)=2a0+n=1[ancos(nx)+bnsin(nx)]

    ②: a n = 2 T ⋅ ∫ − T 2 T 2 f ( x ) ⋅ c o s ( n x ) d x a_n = \frac{2}{T} \cdot \int _{-\frac{T}{2}}^{\frac{T}{2}}f(x)\cdot cos(nx)dx an=T22T2Tf(x)cos(nx)dx

    ③: b n = 2 T ⋅ ∫ − T 2 T 2 f ( x ) ⋅ s i n ( n x ) d x b_n = \frac{2}{T} \cdot \int _{-\frac{T}{2}}^{\frac{T}{2}}f(x)\cdot sin(nx)dx bn=T22T2Tf(x)sin(nx)dx

    顺便提一嘴, a n a_n an求解式中分子的 2 2 2,来源于倍角公式转换引入的,第一项 a 0 a_0 a0没有用这个倍角公式,自然分子没有 2 2 2。(纯属个人瞎掰)

    3.引入复指数

    引入复指数后的推导,在下一章节中详细说明。

    4.总结

    带公式求解,不是天朝的精髓么?可我们为什么要用公式去学习?公式不是为了方便使用才被创造出来的么?理解了,公式才有价值;懂都不懂,公式有什么用?如果仔细看教材,无论教材好坏都有这个的推导。但实际上,56学时的课、5个学分,有几个会真真正正去讲解、去推导?(本段纯属个人发牢骚,我的文章都是这个风格,多包涵)

    三、离散傅里叶变换(DFT)

    上一章介绍了傅里叶级数和傅里叶变换,其实就是利用了一个思想:任何周期函数都可以用正弦函数或者余弦函数的无穷级数来表示。这句话很容易理解,就是证明嘛。用数学公式也是更加简洁:

    f ( t ) = ∑ n = 0 ∞ [ A n ⋅ c o s ( n ω t + ϕ 1 ) ] f(t) = \sum _{n=0}^\infin [A_n \cdot cos(n\omega t + \phi _1) ] f(t)=n=0[Ancos(nωt+ϕ1)]

    或者:

    f ( t ) = ∑ n = 0 ∞ [ B n ⋅ s i n ( n ω t + ϕ 2 ) ] f(t) = \sum _{n=0}^\infin [B_n \cdot sin(n\omega t + \phi _2) ] f(t)=n=0[Bnsin(nωt+ϕ2)]

    可见,我们要求的就是这个 A n A_n An ϕ 1 \phi _1 ϕ1或者 B n B_n Bn ϕ 2 \phi _2 ϕ2

    可能有的人会有疑问了?为什么是 n ω t n\omega t nωt
    其实什么都行,上一章节再推导傅里叶级数时我时从 c o s ( n x ) cos(nx) cos(nx)开始的,然后不断的替换、替换。$n\omega 只 是 表 征 了 只是表征了 cos() 的 频 率 , 也 就 是 周 期 。 你 要 是 高 兴 , 你 用 的频率,也就是周期。你要是高兴,你用 f(t) = \sum _{n=0}^\infin [A_n \cdot cos(狗t + \phi _1) ]$ 都行。
    ω \omega ω,你就可以替换成 ω = 2 π f \omega = 2\pi f ω=2πf;而这个 f f f,不就是我们平常说的多少多少赫兹么。同样,你如果不按照这个替换流程,只不过你最终得到的 f f f可能为 0.01 0.01 0.01赫兹? 2 \sqrt{2} 2 赫兹?这谁顶得住?本身我们得到的就是频谱图,画图而已;况且傅里叶变换的思想就是无限的去逼近,能不能逼近是质的问题,有多逼近是量的问题。而且如果你要是突发奇想用 n E nE nE来表征,说不定你就能用 E = m c 2 E=mc^2 E=mc2来搞事情了。

    扯的有点多,接下来该进入正题了。

    1.离散傅里叶变换的求解

    既然已经知道了傅里叶变换的本质,那么我们不妨就以一个例子来说明:

    f ( t ) = ∑ n = 0 ∞ [ A n ⋅ c o s ( n ω t + ϕ 1 ) ] f(t) = \sum _{n=0}^\infin [A_n \cdot cos(n\omega t + \phi _1) ] f(t)=n=0[Ancos(nωt+ϕ1)]

    这个式子很清楚表示了任何一个周期函数 f ( x ) f(x) f(x)如何被无穷多个余弦函数逼近的。为了求解相对简单,求解的余弦函数 c o s ( n ω x ) cos(n\omega x) cos(nωx)的频率都是整数(如果你愿意,你可以把 ∑ n = 0 ∞ \sum _{n=0}^\infin n=0换成 ∫ n = 0 ∞ \int_{n=0}^\infin n=0,因为积分的本质就是累和,只是累和的步进值无穷小;为什么从0开始?因为没有负频率,当然现在还不去牵扯复指数的那一套东西)

    然后利用我们上一章节推导的,利用三角函数的正交性,积分求系数
    这里我抄过来:

    ①: f ( t ) = a 0 2 + ∑ n = 1 ∞ [ a n ⋅ c o s ( n ω t ) + b n ⋅ s i n ( n ω t ) ] f(t) = \frac{a_0}{2} + \sum _{n=1}^\infin [a_n \cdot cos(n\omega t) + b_n \cdot sin(n\omega t)] f(t)=2a0+n=1[ancos(nωt)+bnsin(nωt)]

    ②: a n = 2 T ⋅ ∫ − T 2 T 2 f ( t ) ⋅ c o s ( n ω t ) d t a_n = \frac{2}{T} \cdot \int _{-\frac{T}{2}}^{\frac{T}{2}}f(t)\cdot cos(n\omega t)dt an=T22T2Tf(t)cos(nωt)dt

    ③: b n = 2 T ⋅ ∫ − T 2 T 2 f ( t ) ⋅ s i n ( n ω t ) d t b_n = \frac{2}{T} \cdot \int _{-\frac{T}{2}}^{\frac{T}{2}}f(t)\cdot sin(n\omega t)dt bn=T22T2Tf(t)sin(nωt)dt

    因为式子①把 A n ⋅ c o s ( n ω t + ϕ 1 ) A_n\cdot cos(n\omega t + \phi _1) Ancos(nωt+ϕ1)分解成了一个 a n ⋅ c o s ( n ω t ) a_n\cdot cos(n\omega t) ancos(nωt)和一个 b n ⋅ s i n ( n ω t ) b_n\cdot sin(n\omega t) bnsin(nωt)

    所以求出来的 a n a_n an b n b_n bn还要转换回去:

    a n ⋅ c o s ( n ω t ) + b n ⋅ s i n ( n ω t ) = a n 2 + b n 2 ⋅ [ a n a n 2 + b n 2 ⋅ c o s ( n ω t ) + b n a n 2 + b n 2 ⋅ s i n ( n ω t ) ] = a n 2 + b n 2 ⋅ c o s ( n ω t + ϕ 1 ) a_n\cdot cos(n\omega t) + b_n\cdot sin(n \omega t) = \sqrt{a_n^2+b_n^2}\cdot [\frac{a_n}{\sqrt{a_n^2+b_n^2}}\cdot cos(n\omega t) + \frac{b_n}{\sqrt{a_n^2+b_n^2}}\cdot sin(n\omega t)] = \sqrt{a_n^2+b_n^2}\cdot cos(n\omega t + \phi_1) ancos(nωt)+bnsin(nωt)=an2+bn2 [an2+bn2 ancos(nωt)+an2+bn2 bnsin(nωt)]=an2+bn2 cos(nωt+ϕ1)

    所以最终的系数为: a n 2 + b n 2 \sqrt{a_n^2+b_n^2} an2+bn2
    相位为: ϕ 1 \phi_1 ϕ1,这个需要大家自己去计算下,我这里就不去计算了

    好了现在问题来了,如何求离散变换下的 a n a_n an b n b_n bn

    下面两个式子是我们已经推导出的连续变换下的 a n a_n an b n b_n bn

    a n = 2 T ⋅ ∫ − T 2 T 2 f ( t ) ⋅ c o s ( n ω t ) d t a_n = \frac{2}{T} \cdot \int _{-\frac{T}{2}}^{\frac{T}{2}}f(t)\cdot cos(n\omega t)dt an=T22T2Tf(t)cos(nωt)dt

    b n = 2 T ⋅ ∫ − T 2 T 2 f ( t ) ⋅ s i n ( n ω t ) d t b_n = \frac{2}{T} \cdot \int _{-\frac{T}{2}}^{\frac{T}{2}}f(t)\cdot sin(n\omega t)dt bn=T22T2Tf(t)sin(nωt)dt

    那么我们不妨设一个周期中采样了16个离散点,即 N = 16 N=16 N=16
    则:

    a n = 2 N ⋅ ∑ 0 N − 1 f ( t N ) ⋅ c o s ( n ω t N ) a_n = \frac{2}{N} \cdot \sum _{0}^{N-1}f(\frac{t}{N})\cdot cos(n\omega \frac{t}{N}) an=N20N1f(Nt)cos(nωNt)

    b n = 2 N ⋅ ∑ 0 N − 1 f ( t N ) ⋅ s i n ( n ω t N ) b_n = \frac{2}{N} \cdot \sum _{0}^{N-1}f(\frac{t}{N})\cdot sin(n\omega \frac{t}{N}) bn=N20N1f(Nt)sin(nωNt)

    要是觉得这个和一些教材上的不一样,我们再变换一下,用 ω = 2 π f \omega = 2\pi f ω=2πf替换:

    a n = 2 N ⋅ ∑ 0 N − 1 f ( t N ) ⋅ c o s ( 2 n π f t N ) a_n = \frac{2}{N} \cdot \sum _{0}^{N-1}f(\frac{t}{N})\cdot cos(2n\pi f\frac{t}{N}) an=N20N1f(Nt)cos(2nπfNt)

    b n = 2 N ⋅ ∑ 0 N − 1 f ( t N ) ⋅ s i n ( 2 n π f t N ) b_n = \frac{2}{N} \cdot \sum _{0}^{N-1}f(\frac{t}{N})\cdot sin(2n\pi f\frac{t}{N}) bn=N20N1f(Nt)sin(2nπfNt)

    然后,然后就没了。

    2.C语言实现

    其实,C语言中实现DFT很简单,就是两层循环:内层循环算这个 ∑ \sum 得到一个 a n a_n an,外层循环算全部 N N N个点的 a n a_n an

    #define PI			3.1415926f
    #define N			8
    
    //@brief:	离散傅里叶变换DFT
    //@param:	*raw_data[in],		原始数据
    //@param:	*dft_coeffi[out],	变换后的幅值序列
    //@param:	len[in],			数据长度
    //@return:	void
    //@note:	输入raw_data和输出dft_result长度相同,都为len
    //@note:	暂未求相位
    //@note:	最高支持N=32768
    void dft(double *raw_data, double *dft_coeffi, int len)
    {
    	double coeffi_cos;
    	double coeffi_sin;
    	
    	for(short freq=0; freq<N; freq++)
    	{
    		coeffi_cos = 0;
    		coeffi_sin = 0;
    		
    		for(short n=0; n<N; n++)
    		{
    			coeffi_cos += raw_data[n] * cos(2*PI*freq*n/N);
    			coeffi_sin += raw_data[n] * sin(2*PI*freq*n/N);
    		}
    		
    		dft_coeffi[freq]  = sqrt(coeffi_cos*coeffi_cos + coeffi_sin*coeffi_sin); //可优化,sin(0)=0
    		dft_coeffi[freq] *= (freq == 0) ? (1.0/N) : (2.0/N);
    	}
    }
    

    但是很多时候我们在网上找到的代码貌似不是这样,很多都是告诉我们:

    1. 先定义一个新的类型:复数类型,就是a+bi或者a+bj
    2. 然后重载一下+*这两个运算符
    3. 然后就巴拉巴拉运算

    实际上,这只用了一些高级语言(如C++)的运算符重载技术,以及引入了复指数、复平面的概念(如果你有留意到,同时这些代码如果写注释的话,一般是这样子的: a n c o s ( 2 π f t N ) − j ⋅ b n ( 2 π f t N ) a_ncos(2\pi f\frac{t}{N}) - j\cdot b_n(2\pi f\frac{t}{N}) ancos(2πfNt)jbn(2πfNt),虚部带个负号)

    但本质上,就是我用C语言中写的那个处理流程,随着 N N N增大,时间复杂度可想而知。但是,这个就是离散傅里叶最原始的样子,暴力求解。时间复杂度很好计算,两层循环, O ( n 2 ) O(n^2) O(n2)。如果还要去区分什么加法、乘法次数,数两次循环的+*个数就行了,但是有一点得提醒一下:你得先理解大O怎么求时间复杂度的;以及实际代码需不需要先被优化,然后再去求

    下面,来解释一下为什么: a n c o s ( 2 π f t N ) − j ⋅ b n ( 2 π f t N ) a_ncos(2\pi f\frac{t}{N}) - j\cdot b_n(2\pi f\frac{t}{N}) ancos(2πfNt)jbn(2πfNt),虚部带个负号。

    3.引入复指数

    最常见的离散傅里叶变换的表达式

    X [ k ] = ∑ n = 0 N − 1 x ( n ) ⋅ e − j 2 π k f n N X[k] = \sum _{n=0}^{N-1} x(n)\cdot e^{-j2\pi kf\frac{n}{N}} X[k]=n=0N1x(n)ej2πkfNn

    如果你有兴趣,直接用欧拉公式 e j x = c o s ( x ) + j s i n ( x ) e^{jx} = cos(x) + jsin(x) ejx=cos(x)+jsin(x)替换展开就可以了

    X [ k ] = ∑ n = 0 N − 1 x ( n ) ⋅ [ c o s ( − 2 π k n N ) + j s i n ( − 2 π k n N ) ] X[k] = \sum _{n=0}^{N-1} x(n)\cdot [cos(-2\pi k \frac{n}{N}) + jsin(-2\pi k \frac{n}{N})] X[k]=n=0N1x(n)[cos(2πkNn)+jsin(2πkNn)]

    利用 c o s ( x ) cos(x) cos(x) s i n ( x ) sin(x) sin(x)的奇偶性

    X [ k ] = ∑ n = 0 N − 1 x ( n ) ⋅ [ c o s ( 2 π k n N ) − j s i n ( 2 π k n N ) ] X[k] = \sum _{n=0}^{N-1} x(n)\cdot [cos(2\pi k \frac{n}{N}) - jsin(2\pi k \frac{n}{N})] X[k]=n=0N1x(n)[cos(2πkNn)jsin(2πkNn)]

    然后再告诉你,这个 j j j是一个表示旋转的标识,和 − 1 -1 1一样。在复平面里, j j j标识逆时针转 90 ° 90° 90° − 1 -1 1就表示逆时针转 180 ° 180° 180°,因为 − 1 = j × j -1 = j \times j 1=j×j。总之 j j j在表达式中的存在只是为了告诉我们:我后面的数据是被旋转了的,你要注意。所以我们在处理时没必要管这个 j j j

    再因为最后我们求的 X [ k ] = a k 2 + b k 2 X[k]=\sqrt{a_k^2 + b_k^2} X[k]=ak2+bk2 ,平方和求开根号,所以 X [ k ] = a k 2 + ( − b k ) 2 = a k 2 + b k 2 X[k]=\sqrt{a_k^2 + (-b_k)^2} = \sqrt{a_k^2 + b_k^2} X[k]=ak2+(bk)2 =ak2+bk2

    所以,可以得出一个结论:DFT在引入复指数、复平面概念后去求 c o s ( x ) cos(x) cos(x) s i n ( x ) sin(x) sin(x)的系数,所以用绝对值就可以了,系数的正负是无所谓的。

    4.什么是复平面

    关于复平面,我直接引用一位大神slowmovingsnail的CSDN文章

    傅里叶分析的方方面面:复正弦、负频率

    这张图,如果幸运的话,我可能还会讲到,不过现在只是用它来说明一个问题:通过欧拉公式的替换,把 c o s ( x ) cos(x) cos(x)或者 s i n ( x ) sin(x) sin(x)用两个复数来表示

    即:

    c o s ( θ ) = 1 2 ⋅ [ e j θ + e − j θ ] cos(\theta) = \frac{1}{2} \cdot [e^{j\theta} + e^{-j\theta}] cos(θ)=21[ejθ+ejθ]

    s i n ( θ ) = 1 2 j ⋅ [ e j θ − e − j θ ] sin(\theta) = \frac{1}{2j} \cdot [e^{j\theta} - e^{-j\theta}] sin(θ)=2j1[ejθejθ]

    而我们已经知道,任何周期函数都可以展开成:
    这里我就不单独写 a 0 a_0 a0了,反正既然能看到这,想必也都知道为什么了

    f ( t ) = ∫ n = 0 ∞ [ a n ⋅ c o s ( n ω t ) + b n ⋅ s i n ( n ω t ) ] f(t) = \int _{n=0}^\infin [a_n \cdot cos(n\omega t) + b_n \cdot sin(n\omega t)] f(t)=n=0[ancos(nωt)+bnsin(nωt)]

    带入替换:

    f ( t ) = ∫ n = 0 ∞ [ a n ⋅ 1 2 ⋅ ( e j n ω t + e − j n ω t ) + b n ⋅ 1 2 j ⋅ ( e j n ω t − e − j n ω t ) ] f(t) = \int _{n=0}^\infin [a_n \cdot \frac{1}{2} \cdot (e^{jn\omega t} + e^{-jn\omega t}) + b_n \cdot \frac{1}{2j} \cdot (e^{jn\omega t} - e^{-jn\omega t})] f(t)=n=0[an21(ejnωt+ejnωt)+bn2j1(ejnωtejnωt)]

    整理:

    f ( t ) = ∫ n = 0 ∞ [ a n ⋅ 1 2 ⋅ ( e j n ω t + e − j n ω t ) + b n ⋅ 1 2 j ⋅ ( e j n ω t − e − j n ω t ) ] f(t) = \int _{n=0}^\infin [a_n \cdot \frac{1}{2} \cdot (e^{jn\omega t} + e^{-jn\omega t}) + b_n \cdot \frac{1}{2j} \cdot (e^{jn\omega t} - e^{-jn\omega t})] f(t)=n=0[an21(ejnωt+ejnωt)+bn2j1(ejnωtejnωt)]

    说实话,我也很难看出还能怎么往下做。但我只能说我们的目的是尽可能转化成只有一个 e j θ e^{j\theta} ejθ的形式,所以,凑就完了。

    f ( t ) = ∫ n = 0 ∞ [ a n 2 ⋅ ( e j n ω t + e − j n ω t ) + b n 2 ⋅ ( − j ) ⋅ ( e j n ω t − e − j n ω t ) ] f(t) = \int _{n=0}^\infin [\frac{a_n}{2} \cdot (e^{jn\omega t} + e^{-jn\omega t}) + \frac{b_n }{2}\cdot (-j) \cdot (e^{jn\omega t} - e^{-jn\omega t})] f(t)=n=0[2an(ejnωt+ejnωt)+2bn(j)(ejnωtejnωt)]

    上式的处理中,利用了: 1 = j 4 = j 2 ⋅ j ⋅ j = ( − j ) ⋅ j 1 = j^4 = j^2 \cdot j \cdot j = (-j) \cdot j 1=j4=j2jj=(j)j

    所以,再整理:

    f ( t ) = ∫ n = 0 ∞ [ ( a n − j b n ) 2 ⋅ e j n ω t + ( a n + j b n ) 2 ⋅ e − j n ω t ] f(t) = \int _{n=0}^\infin [\frac{(a_n - jb_n)}{2} \cdot e^{jn\omega t} + \frac{(a_n + jb_n)}{2} \cdot e^{-jn\omega t}] f(t)=n=0[2(anjbn)ejnωt+2(an+jbn)ejnωt]

    令: c n 1 = ( a n − j b n ) 2 c_{n1} = \frac{(a_n - jb_n)}{2} cn1=2(anjbn) c n 2 = ( a n + j b n ) 2 c_{n2} = \frac{(a_n + jb_n)}{2} cn2=2(an+jbn)

    可以发现, c n 1 c_{n1} cn1 c n 2 c_{n2} cn2是共轭的(即模相等,虚部相反)

    f ( t ) = ∫ n = 0 ∞ [ c n 1 ⋅ e j n ω t + c n 2 ⋅ e − j n ω t ] = ∑ − ∞ + ∞ ∣ d n ∣ ⋅ e j n ω t f(t) = \int _{n=0}^\infin [c_{n1} \cdot e^{jn\omega t} + c_{n2} \cdot e^{-jn\omega t}] = \sum _{-\infin}^{+\infin} |d_n| \cdot e^{jn\omega t} f(t)=n=0[cn1ejnωt+cn2ejnωt]=+dnejnωt注意,这里先这么写,这里不是最终的!

    上一步的推导好像理解起来没问题,可又好像有问题。上一步还是 a n − j b n a_n - jb_n anjbn,下一步直接 ∣ d n ∣ |d_n| dn了, j j j直接没了?

    其实,这里借助了一个概念:旋转因子
    旋转因子就是 e j n ω t e^{jn\omega t} ejnωt e − j n ω t e^{-jn\omega t} ejnωt,前者表示逆时针转 n ω t n\omega t nωt,后者表示顺时针转 n ω t n\omega t nωt

    再加上: c n 1 = ( a n − j b n ) 2 c_{n1} = \frac{(a_n - jb_n)}{2} cn1=2(anjbn) c n 2 = ( a n + j b n ) 2 c_{n2} = \frac{(a_n + jb_n)}{2} cn2=2(an+jbn)是共轭的,也就是在复平面内, c n 1 c_{n1} cn1 c n 2 c_{n2} cn2这两个向量关于横轴(也就是实轴)对称(我们这里称 c n 1 c_{n1} cn1 c n 2 c_{n2} cn2为复平面里的两个向量,绕原点旋转的向量)。再加上旋转因子 e j n ω t e^{jn\omega t} ejnωt e − j n ω t e^{-jn\omega t} ejnωt分别表示向相反方向旋转相同的角度,所以 c n 1 ⋅ e j n ω t c_{n1} \cdot e^{jn\omega t} cn1ejnωt c n 2 ⋅ e − j n ω t c_{n2} \cdot e^{-jn\omega t} cn2ejnωt也是关于横轴(实轴)对称。
    然后我们把 c n 1 ⋅ e j n ω t c_{n1} \cdot e^{jn\omega t} cn1ejnωt c n 2 ⋅ e − j n ω t c_{n2} \cdot e^{-jn\omega t} cn2ejnωt这两个向量的模长拿出来,就是 ∣ d n ∣ |d_n| dn,其中 ∣ d n ∣ = ∣ c n 1 ∣ = ∣ c n 2 ∣ |d_n| = |c_{n1}| = |c_{n2}| dn=cn1=cn2。但是这样就丢失了 c n 1 c_{n1} cn1 c n 2 c_{n2} cn2本身的角度信息,所以我们就需要把丢失的角度补到后面,即: ∣ d n ∣ e j ( n ω t + θ c n ) |d_n|e^{j(n\omega t + \theta _{c_n})} dnej(nωt+θcn) ∣ d n ∣ e − j ( n ω t + θ c n ) |d_n|e^{-j(n\omega t + \theta _{c_n})} dnej(nωt+θcn)

    所以,我们把前面挖的坑填上:

    f ( t ) = ∫ n = 0 ∞ [ c n 1 ⋅ e j n ω t + c n 2 ⋅ e − j n ω t ] = ∑ − ∞ + ∞ ∣ d n ∣ ⋅ e j ( n ω t + θ c n ) f(t) = \int _{n=0}^\infin [c_{n1} \cdot e^{jn\omega t} + c_{n2} \cdot e^{-jn\omega t}] = \sum _{-\infin}^{+\infin} |d_n| \cdot e^{j(n\omega t + \theta _{c_n})} f(t)=n=0[cn1ejnωt+cn2ejnωt]=+dnej(nωt+θcn)

    当然,也可以从图形结合来理解:

    1 2 a n + j 1 2 b n = > ∣ c n ∣ e j θ c n \frac{1}{2}a_n + j\frac{1}{2}b_n => |c_n|e^{j\theta _{c_n}} 21an+j21bn=>cnejθcn,注意我这里用的不是一个等于号,我只是想表达一个意思: 1 2 a n + j 1 2 b n \frac{1}{2}a_n + j\frac{1}{2}b_n 21an+j21bn可以在复平面中找到一个绕原点旋转 θ c n \theta _{c_n} θcn的向量来对应,这个向量的长度为 ∣ c n ∣ |c_n| cn借助图形可以很容易分析出来,只是我写这些公式太麻烦了。

    好了,通过以上的分析,我们就得到了引入复指数、复平面后的展开式:

    f ( t ) = ∫ − ∞ + ∞ ∣ d n ∣ ⋅ e j ( n ω t + θ c n ) f(t) = \int _{-\infin}^{+\infin} |d_n| \cdot e^{j(n\omega t + \theta _{c_n})} f(t)=+dnej(nωt+θcn)

    但更多的时候,我们关注的是这个幅值: ∣ d n ∣ |d_n| dn
    同样,上式是连续函数的展开,离散函数的展开如下

    f [ n ] = ∑ − ∞ + ∞ ∣ d n ∣ ⋅ e j ( k ω n N + θ c n ) f[n] = \sum _{-\infin}^{+\infin} |d_n| \cdot e^{j(k\omega \frac{n}{N} + \theta _{c_n})} f[n]=+dnej(kωNn+θcn)

    但是,展开式中带着 θ c n \theta _{c_n} θcn,在求解时很麻烦。

    既然这样,我们不妨从复指数展开式倒推回三角展开式

    只不过,用的: f [ n ] = ∑ − ∞ + ∞ ∣ d n ∣ ⋅ e j ( k ω n N ) f[n] = \sum _{-\infin}^{+\infin} |d_n| \cdot e^{j(k\omega \frac{n}{N})} f[n]=+dnej(kωNn)

    相当于,我们只用了 c o s ( x ) cos(x) cos(x)来展开

    即: f [ n ] = ∑ k = 0 ∞ a n ⋅ c o s ( k ω n N ) = ∑ k = 0 ∞ a n 2 ⋅ [ e j k ω n N + e − j k ω n N ] = ∑ − ∞ + ∞ ∣ c n ∣ ⋅ e j k ω n N f[n] = \sum _{k=0}^{\infin} a_n \cdot cos(k\omega \frac{n}{N}) = \sum _{k=0}^{\infin} \frac{a_n}{2}\cdot [e^{jk\omega \frac{n}{N}} + e^{-jk\omega \frac{n}{N}}] = \sum _{-\infin}^{+\infin} |c_n| \cdot e^{jk\omega \frac{n}{N}} f[n]=k=0ancos(kωNn)=k=02an[ejkωNn+ejkωNn]=+cnejkωNn

    这样就是说:我们假定是用 c o s ( k ω ) cos(k \omega) cos(kω)来展开,而且我们丢失拟合精度,因为正统的傅里叶级数是由两部分构成的, c o s ( k ω ) cos(k \omega) cos(kω) s i n ( k ω ) sin(k \omega) sin(kω),我们这里丢失了 s i n ( k ω ) sin(k \omega) sin(kω)部分。实际上,我实在无法用这个说法来说服自己,可能是我推导错了

    我目前无法解决的问题是:为什么从三角推导过来的复指数展开式 f [ n ] = ∑ − ∞ + ∞ ∣ d n ∣ ⋅ e j ( k ω n N + θ c n ) f[n] = \sum _{-\infin}^{+\infin} |d_n| \cdot e^{j(k\omega \frac{n}{N} + \theta _{c_n})} f[n]=+dnej(kωNn+θcn) 明明有这个 θ c n \theta _{c_n} θcn,可实际的展开式却没有 f [ n ] = ∑ − ∞ + ∞ ∣ d n ∣ ⋅ e j ( k ω n N ) f[n] = \sum _{-\infin}^{+\infin} |d_n| \cdot e^{j(k\omega \frac{n}{N})} f[n]=+dnej(kωNn)?如果有答案,还请不吝赐教,谢谢。

    5.为什么DFT表达式里是 − j 2 π k f n N -j2\pi kf\frac{n}{N} j2πkfNn

    把DFT表达式抄过来镇楼:

    X [ k ] = ∑ n = 0 N − 1 x ( n ) ⋅ e − j 2 π k f n N X[k] = \sum _{n=0}^{N-1} x(n)\cdot e^{-j2\pi kf\frac{n}{N}} X[k]=n=0N1x(n)ej2πkfNn

    再抄一个离散傅里叶展开式来镇二楼:

    x [ n ] = ∑ − ∞ + ∞ X [ k ] ⋅ e j 2 π k f n N x[n] = \sum _{-\infin}^{+\infin} X[k] \cdot e^{j2\pi kf \frac{n}{N}} x[n]=+X[k]ej2πkfNn

    冷不丁看这两个式子很像,但实际上差的很远,不仅仅是一个是DFT,一个是展开,主要是两个式子的 ∑ \sum 是不同的:

    1. DFT中 ∑ \sum 是指对一个周期的所有采样点进行处理,即步进变量是 n n n
    2. 展开式中 ∑ \sum 是指的所有频率,即步进变量是 k k k

    利用正交性,再对展开式两边的 n n n积分累和,就可以得到DFT:

    ∑ n = 0 N − 1 x [ n ] ⋅ e j 2 π k f n N = ∑ n = 0 N − 1 { ∑ − ∞ + ∞ X [ k ] ⋅ e j 2 π k f n N } ⋅ e j 2 π k f n N \sum _{n=0}^{N-1} x[n]\cdot e^{j2\pi kf \frac{n}{N}} = \sum _{n=0}^{N-1} \{\sum _{-\infin}^{+\infin} X[k] \cdot e^{j2\pi kf \frac{n}{N}}\}\cdot e^{j2\pi kf \frac{n}{N}} n=0N1x[n]ej2πkfNn=n=0N1{+X[k]ej2πkfNn}ej2πkfNn

    X [ k ] X[k] X[k]是一个待求的值,不包含被积量,所以提出:

    ∑ n = 0 N − 1 x [ n ] ⋅ e j 2 π k f n N = X [ k ] ⋅ ∑ n = 0 N − 1 { ∑ − ∞ + ∞ e j 2 π k f n N } ⋅ e j 2 π k f n N \sum _{n=0}^{N-1} x[n]\cdot e^{j2\pi kf \frac{n}{N}} = X[k] \cdot \sum _{n=0}^{N-1} \{\sum _{-\infin}^{+\infin} e^{j2\pi kf \frac{n}{N}}\}\cdot e^{j2\pi kf \frac{n}{N}} n=0N1x[n]ej2πkfNn=X[k]n=0N1{+ej2πkfNn}ej2πkfNn

    ∑ n = 0 N − 1 x [ n ] ⋅ e j 2 π k f n N = X [ k ] ⋅ e j 2 π k f n N ⋅ e j 2 π k f n N \sum _{n=0}^{N-1} x[n]\cdot e^{j2\pi kf \frac{n}{N}} = X[k] \cdot e^{j2\pi kf \frac{n}{N}} \cdot e^{j2\pi kf \frac{n}{N}} n=0N1x[n]ej2πkfNn=X[k]ej2πkfNnej2πkfNn

    整理,就得到:

    X [ k ] = ∑ n = 0 N − 1 x [ n ] ⋅ e − j 2 π k f n N X[k] = \sum _{n=0}^{N-1} x[n]\cdot e^{-j2\pi kf \frac{n}{N}} X[k]=n=0N1x[n]ej2πkfNn

    所以,好像这一小节我推导了个寂寞。

    简单用复指数证明一下正交性:

    用欧拉公式替换成三角表示,用三角函数的正交性来解释。偷懒ing

    6.总结

    其实这里总结的不仅仅是离散傅里叶变换,实际上我把复指数表示方式的证明也放在了本章节。

    1. 引入复指数、复平面概念后,引入了负频率,这个在实际生活中是没有意义的,因为你找不到频率<0(或者说,在目前的认知中,尚且无法解释负频率的实际物理意义,只能暂时认为是数学运算产生的中间过程量)
    2. 很无奈,有一个问题至今无法解释,实在是个人能力有限

    四、快速傅里叶变换(FFT)

    未完待续。

    五、总结

    世界上的每一个人,都被同一句谎话欺骗着:站在巨人的肩膀上看世界,会看的更远。我们都牢牢记住了这句话,确实看的更远,不过也仅此而已。

    最后,再引用知乎leinlin文章的最后的一句话作为结束:

    我们眼中的世界就像皮影戏的大幕布,幕布的后面有无数的齿轮,大齿轮带动小齿轮,小齿轮再带动更小的。在最外面的小齿轮上有一个小人——那就是我们自己。我们只看到这个小人毫无规律的在幕布前表演,却无法预测他下一步会去哪。而幕布后面的齿轮却永远一直那样不停的旋转,永不停歇。这样说来有些宿命论的感觉。说实话,这种对人生的描绘是我一个朋友在我们都是高中生的时候感叹的,当时想想似懂非懂,直到有一天我学到了傅里叶级数……

    其实,懂得人都懂了,不懂的还是不懂。

    展开全文
  • 1.用Matlab符号运算求解法求单边指数信号f(t)=e−2tu(t)f(t)=e^{-2t}u(t)f(t)=e−2tu(t)的FT MATLAB源程序为: ft=sym('exp(-2*t)*heaviside(t)'); fw=fourier(ft) 运行结果为: fw = 1/(2 + w*1i) 2.用Matlab符号...

    1.用Matlab符号运算求解法求单边指数信号 f ( t ) = e − 2 t u ( t ) f(t)=e^{-2t}u(t) f(t)=e2tu(t)的FT

    MATLAB源程序为:

    ft=sym('exp(-2*t)*heaviside(t)');
    fw=fourier(ft)
    

    运行结果为:

    fw =
    1/(2 + w*1i)
    

    2.用Matlab符号运算求解法求 F ( j w ) = 1 1 + w 2 F(jw)=\frac{1}{1+w^2} F(jw)=1+w21的IFT。

    解:matlab源程序为

    syms t
    fw=sym('1/(1+w^2)');
    ft=ifourier(fw,t)
    

    运行结果为:

    ft =
    exp(-abs(t))/2
    

    3.用Matlab命令绘出信号 f ( t ) = e − 2 t u ( t ) f(t)=e^{-2t}u(t) f(t)=e2tu(t)的频谱图。

    Matlab源程序为:

    ft=sym('exp(-2*t)*heaviside(t)');
    fw=fourier(ft);
    subplot(211)
    ezplot(abs(fw)),grid on
    title('幅度谱')
    phase=atan(imag(fw)/real(fw));
    subplot(212)
    ezplot(phase);grid on
    title('相位谱')
    

    调用heaviside函数需要与我们定义的函数名一致,一开始h大写了,就报错,后来改为小写就正确了。
    在这里插入图片描述

    展开全文
  • 傅立叶变换MATLAB的实现及傅里叶变换性质的分析。实验目的: 1.利用MATLAB分析非周期信号的频谱 2.观察信号频谱变化验证傅里叶变换性质

    【 实验目的】

    1.利用MATLAB分析非周期信号的频谱
    2.观察信号频谱变化验证傅里叶变换性质

    【 实验内容】

    在这里插入图片描述

    【 实验报告要求】

    (1)记录实验一和实验三中的波形;
    (2)总结实验二中频谱特性曲线变化的特点;
    (3)实验目的和实验过程进行总结。

    Matlab程序一:

    f=sym('exp(-2*t)*heaviside(t)');%方程
    F=fourier(f);%傅里叶变换
    FT=subs(F,[-2*pi:0.01*pi:2*pi]);%subs()函数来得到傅立叶变换的数值解
    subplot(4,2,1); %作图区域划分    
    plot([-2*pi:0.01*pi:2*pi],abs(FT));%绘制幅频特性曲线
    title('幅频曲线-Make by 磊');
    xlabel('w');ylabel('幅度');
    subplot(4,2,2);%作图区域划分  
    plot([-2*pi:0.01*pi:2*pi],angle(FT));%绘制相频特性曲线
    title('相频曲线-Make by 磊');
    xlabel('w');ylabel('相位(弧度)');grid on; 
    f=sym('exp(-3*abs(t))'); %方程
    F=fourier(f);%傅里叶变换
    FT=subs(F,[-2*pi:0.01*pi:2*pi]);%subs()函数来得到傅立叶变换的数值解
    subplot(4,2,3);%作图区域划分     
    plot([-2*pi:0.01*pi:2*pi],abs(FT));%绘制幅频特性曲线
    title('幅频曲线-Make by 磊');
    xlabel('w');ylabel('幅度');
    subplot(4,2,4);%作图区域划分   
    plot([-2*pi:0.01*pi:2*pi],angle(FT));%绘制相频特性曲线
    title('相频曲线-Make by 磊');
    xlabel('w');ylabel('相位(弧度)');grid on;
    f=sym('heaviside(t+2)-heaviside(t-2)');%方程 
    F=fourier(f);%傅里叶变换
    t1=-10*pi:0.01*pi:-0.001;
    t2=0.001*pi:0.01*pi:10*pi;
    t=[t1 t2];
    FT=subs(F,t);%subs()函数来得到傅立叶变换的数值解
    subplot(4,2,5);%作图区域划分    
    plot(t,abs(FT));%绘制幅频特性曲线
    title('幅频曲线-Make by 磊');
    xlabel('w');ylabel('幅度');grid on;
    subplot(4,2,6);%作图区域划分   
    plot(t,angle(FT));%绘制相频特性曲线
    title('相频曲线-Make by 磊');
    xlabel('w');ylabel('相位(弧度)');grid on;
    f=sym('sin(2*t)/(2*pi*t)');%方程 
    F=fourier(f);%傅里叶变换
    FT=subs(F,[-2*pi:0.01*pi:2*pi]);%subs()函数来得到傅立叶变换的数值解
    subplot(4,2,7);%作图区域划分     
    plot([-2*pi:0.01*pi:2*pi],abs(FT));%绘制幅频特性曲线
    title('幅频曲线-Make by 磊');
    xlabel('w');ylabel('幅度');
    subplot(4,2,8);%作图区域划分   
    plot([-2*pi:0.01*pi:2*pi],angle(FT));%绘制相频特性曲线
    title('相频曲线-Make by 磊');
    xlabel('w');ylabel('相位(弧度)');grid on;
    
    

    运行结果一:

    在这里插入图片描述

    Matlab程序二:

    f=sym('exp(-2*t)*heaviside(t)');%方程
    F=fourier(f);%傅里叶变换
    FT=subs(F,[-2*pi:0.01*pi:2*pi]);%subs()函数来得到傅立叶变换的数值解
    subplot(3,2,1);%作图区域划分     
    plot([-2*pi:0.01*pi:2*pi],abs(FT));%绘制幅频特性曲线
    title('幅频曲线-Make by 磊');
    xlabel('w');ylabel('幅度');
    subplot(3,2,2); %作图区域划分  
    plot([-2*pi:0.01*pi:2*pi],angle(FT));%绘制相频特性曲线
    title('相频曲线-Make by 磊');
    xlabel('w');ylabel('相位(弧度)');grid on;
    f=sym('exp(-2*(t-0.6))*heaviside(t-0.6)'); %方程
    F=fourier(f);%傅里叶变换
    FT=subs(F,[-2*pi:0.01*pi:2*pi]);%subs()函数来得到傅立叶变换的数值解
    subplot(3,2,3);%作图区域划分     
    plot([-2*pi:0.01*pi:2*pi],abs(FT));%绘制幅频特性曲线
    title('幅频曲线-Make by 磊');
    xlabel('w');ylabel('幅度');
    subplot(3,2,4);%作图区域划分   
    plot([-2*pi:0.01*pi:2*pi],angle(FT));%绘制相频特性曲线
    title('相频曲线-Make by 磊');
    xlabel('w');ylabel('相位(弧度)');grid on;
    f=sym('exp(-2*t)*heaviside(t)*exp(i*2*t)');%方程
    F=fourier(f);%傅里叶变换
    FT2=subs(F,[-2*pi:0.01*pi:2*pi]);%subs()函数来得到傅立叶变换的数值解
    subplot(3,2,5);%作图区域划分
    plot([-2*pi:0.01*pi:2*pi],abs(FT2));%绘制幅频特性曲线
    grid on;
    title('幅频曲线-Make by 磊');
    xlabel('w');ylabel('幅度');
    subplot(3,2,6);%作图区域划分     plot([-2*pi:0.01*pi:2*pi],angle(FT2));%绘制相频特性曲线
    grid on;
    title('相频曲线-Make by 磊');
    xlabel('w');ylabel('相位(弧度)');
    
    

    运行结果二:

    在这里插入图片描述

    Matlab程序三:

    f=sym('sin(2*t)/(2*pi*t)');
    F=fourier(f);%傅里叶变换
    FT=subs(F,[-2*pi:0.01*pi:2*pi]);%subs()函数来得到傅立叶变换的数值解
    subplot(2,2,1);%作图区域划分    
    plot([-2*pi:0.01*pi:2*pi],abs(FT));%绘制幅频特性曲线
    title('幅频曲线-Make by 磊');
    xlabel('w');ylabel('幅度');
    subplot(2,2,2);%作图区域划分  
    plot([-2*pi:0.01*pi:2*pi],angle(FT));%绘制相频特性曲线
    title('相频曲线-Make by 磊');
    xlabel('w');ylabel('相位(弧度)');grid on;
    f=sym('(sin(2*t)/(2*pi*t))*cos(6*t)'); %方程
    F=fourier(f);%傅里叶变换
    FT=subs(F,[-2*pi:0.01*pi:2*pi]);%subs()函数来得到傅立叶变换的数值解
    subplot(2,2,3);%作图区域划分     
    plot([-2*pi:0.01*pi:2*pi],abs(FT));%绘制幅频特性曲线
    title('幅频曲线-Make by 磊');
    xlabel('w');ylabel('幅度');
    subplot(2,2,4);%作图区域划分  
    plot([-2*pi:0.01*pi:2*pi],angle(FT));%绘制相频特性曲线
    title('相频曲线-Make by 磊');
    xlabel('w');ylabel('相位(弧度)');grid on;
    
    

    运行结果三:

    在这里插入图片描述

    Matlab程序四:

    f=sym('1/pi*(heaviside(t+0.5*pi)-heaviside(t-0.5*pi))');%方程 
    F=fourier(f);%傅里叶变换
    t1=-10*pi:0.01*pi:-0.001;