从高等数学课本中我们了解到如果一个周期为T的周期函数 f ( t ) f(t) f(t)满足收敛定理的条件，则这个函数可以被展开为傅里叶级数：
- f ( x ) = a 0 + ∑ n = 1 ∞ [ a n c o s ( n w 1 t ) + b n s i n ( n w 1 t ) ] f(x)=a_0+\sum\limits_{n=1}^{\infty}[a_ncos(nw_1t)+b_nsin(nw_1t)] f(x)=a0+n=1∑∞[ancos(nw1t)+bnsin(nw1t)]， w 1 = 2 π T w_1=\frac{2\pi} {T} w1=T2π为角频率，等式一。
- a 0 = 1 T ∫ t 0 t 0 + T f ( t ) d t a_0=\frac{1}{T}\int_{t_0}^{t_0+T}f(t)d_t a0=T1∫t0t0+Tf(t)dt为直流分量即频率为0的分量。
- a n = 2 T ∫ t 0 t 0 + T f ( t ) c o s ( n w 1 t ) d t a_n=\frac{2}{T}\int_{t_0}^{t_0+T}f(t)cos(nw_1t)d_t an=T2∫t0t0+Tf(t)cos(nw1t)dt为余弦分量的幅度， n = 1 , 2 , . . . , n=1,2,..., n=1,2,...,。
- b n = 2 T ∫ t 0 t 0 + T f ( t ) s i n ( n w 1 t ) d t b_n=\frac{2}{T}\int_{t_0}^{t_0+T}f(t)sin(nw_1t)d_t bn=T2∫t0t0+Tf(t)sin(nw1t)dt为正弦分量的幅度， n = 1 , 2 , . . . , n=1,2,..., n=1,2,...,。
又因为 a n c o s ( n w 1 t ) + b n s i n ( n w 1 t ) = a n 2 + b n 2 [ a n a n 2 + b n 2 c o s ( n w 1 t ) + b n a n 2 + b n 2 s i n ( n w 1 t ) ] a_ncos(nw_1t)+b_nsin(nw_1t)=\sqrt{a_n^2+b_n^2}[\frac{a_n}{\sqrt{a_n^2+b_n^2}}cos(nw_1t)+\frac{b_n}{\sqrt{a_n^2+b_n^2}}sin(nw_1t)] ancos(nw1t)+bnsin(nw1t)=an2+bn2[an2+bn2ancos(nw1t)+an2+bn2bnsin(nw1t)]= a n 2 + b n 2 [ a n a n 2 + b n 2 c o s ( n w 1 t ) − − b n a n 2 + b n 2 s i n ( n w 1 t ) ] \sqrt{a_n^2+b_n^2}[\frac{a_n}{\sqrt{a_n^2+b_n^2}}cos(nw_1t)-\frac{-b_n}{\sqrt{a_n^2+b_n^2}}sin(nw_1t)] an2+bn2[an2+bn2ancos(nw1t)−an2+bn2−bnsin(nw1t)]，从以上等式我们可以推导出：
- f ( t ) = a 0 + a n 2 + b n 2 ∑ n = 1 ∞ s i n ( n w 1 t + θ n ) f(t)=a_0+\sqrt{a_n^2+b_n^2}\sum\limits_{n=1}^{\infty}sin(nw_1t+\theta_n) f(t)=a0+an2+bn2n=1∑∞sin(nw1t+θn)， s i n ( θ n ) = a n a n 2 + b n 2 sin(\theta_n)=\frac{a_n}{\sqrt{a_n^2+b_n^2}} sin(θn)=an2+bn2an， c o s ( θ n ) = b n a n 2 + b n 2 cos(\theta_n)=\frac{b_n}{\sqrt{a_n^2+b_n^2}} cos(θn)=an2+bn2bn
- f ( t ) = a 0 + a n 2 + b n 2 ∑ n = 1 ∞ c o s ( n w 1 t + φ n ) f(t)=a_0+\sqrt{a_n^2+b_n^2}\sum\limits_{n=1}^{\infty}cos(nw_1t+\varphi_n) f(t)=a0+an2+bn2n=1∑∞cos(nw1t+φn)， s i n ( φ n ) = − b n a n 2 + b n 2 sin(\varphi_n)=\frac{-b_n}{\sqrt{a_n^2+b_n^2}} sin(φn)=an2+bn2−bn， c o s ( φ n ) = a n a n 2 + b n 2 cos(\varphi_n)=\frac{a_n}{\sqrt{a_n^2+b_n^2}} cos(φn)=an2+bn2an
从以上的推导中我们可以发现信号 f ( x ) f(x) f(x)被分解为了一系列频率不同( w = n w 1 ， n = 0 , 1 , 2 , . . . , w=nw_1，n=0,1,2,..., w=nw1，n=0,1,2,...,)、幅度( a n 2 + b n 2 ， n = 0 , 1 , 2 , . . . , \sqrt{a_n^2+b_n^2}，n=0,1,2,..., an2+bn2，n=0,1,2,...,)不同和相位( φ n 或 θ n ， n = 0 , 1 , 2 , . . . , \varphi_n或\theta_n，n=0,1,2,..., φn或θn，n=0,1,2,...,)不同的正弦或余弦函数。这就相当于是傅里叶分析的初级形式了，我们这时从时间域转换到了频率域。我们此时可以以频率( w = n w 1 ， n = 0 , 1 , 2 , . . . , w=nw_1，n=0,1,2,..., w=nw1，n=0,1,2,...,)为横坐标，幅度( a n 2 + b n 2 ， n = 0 , 1 , 2 , . . . , \sqrt{a_n^2+b_n^2}，n=0,1,2,..., an2+bn2，n=0,1,2,...,)和相位( φ n 或 θ n ， n = 0 , 1 , 2 , . . . , \varphi_n或\theta_n，n=0,1,2,..., φn或θn，n=0,1,2,...,)为纵坐标作图，这样就更显示了频率域的意义了。一个简单的示意图如图1所示（注意这里的图示不是上面的结果，只是一个简单的示意而已）。
欧拉公式：
- e i x = c o s ( x ) + i ∗ s i n ( x ) e^{ix}=cos(x)+i*sin(x) eix=cos(x)+i∗sin(x)， e − i x = c o s ( x ) − i ∗ s i n ( x ) e^{-ix}=cos(x)-i*sin(x) e−ix=cos(x)−i∗sin(x) (1)
- c o s ( x ) = e i x + e − i x 2 cos(x)=\frac{e^{ix}+e^{-ix}}{2} cos(x)=2eix+e−ix， s i n ( x ) = e i x − e − i x 2 i sin(x)=\frac{e^{ix}-e^{-ix}}{2i} sin(x)=2ieix−e−ix (2)
根据欧拉等式我们可以得到：
- s i n ( n w 1 t + θ n ) sin(nw_1t+\theta_n) sin(nw1t+θn)= e i ( n w 1 t + θ n ) − e − i ( n w 1 t + θ n ) 2 i \frac{e^{i(nw_1t+\theta_n)}-e^{-i(nw_1t+\theta_n)}}{2i} 2iei(nw1t+θn)−e−i(nw1t+θn) = e i n w 1 t e i θ n − e − i n w 1 t e − i θ n 2 i \frac{e^{inw_1t}e^{i\theta_n}-e^{-inw_1t}e^{-i\theta_n}}{2i} 2ieinw1teiθn−e−inw1te−iθn= − i e i n w 1 t e i θ n + i e − i n w 1 t e − i θ n 2 \frac{-ie^{inw_1t}e^{i\theta_n}+ie^{-inw_1t}e^{-i\theta_n}}{2} 2−ieinw1teiθn+ie−inw1te−iθn
- c o s ( n w 1 t + φ n ) cos(nw_1t+\varphi_n) cos(nw1t+φn)= e i ( n w 1 t + φ n ) + e − i ( n w 1 t + φ n ) 2 \frac{e^{i(nw_1t+\varphi_n)}+e^{-i(nw_1t+\varphi_n)}}{2} 2ei(nw1t+φn)+e−i(nw1t+φn) = e i n w 1 t e i φ n + e − i n w 1 t e − i φ n 2 \frac{e^{inw_1t}e^{i\varphi_n}+e^{-inw_1t}e^{-i\varphi_n}}{2} 2einw1teiφn+e−inw1te−iφn
又因为：
- e i θ n = c o s ( θ n ) + i ∗ s i n ( θ n ) e^{i\theta_n}=cos(\theta_n)+i*sin(\theta_n) eiθn=cos(θn)+i∗sin(θn)= b n a n 2 + b n 2 \frac{b_n}{\sqrt{a_n^2+b_n^2}} an2+bn2bn + i a n a n 2 + b n 2 +i\frac{a_n}{\sqrt{a_n^2+b_n^2}} +ian2+bn2an， e − i θ n = c o s ( θ n ) − i ∗ s i n ( θ n ) e^{-i\theta_n}=cos(\theta_n)-i*sin(\theta_n) e−iθn=cos(θn)−i∗sin(θn)= b n a n 2 + b n 2 \frac{b_n}{\sqrt{a_n^2+b_n^2}} an2+bn2bn − i a n a n 2 + b n 2 -i\frac{a_n}{\sqrt{a_n^2+b_n^2}} −ian2+bn2an，
- e i φ n = c o s ( φ n ) + i ∗ s i n ( φ n ) e^{i\varphi_n}=cos(\varphi_n)+i*sin(\varphi_n) eiφn=cos(φn)+i∗sin(φn)= a n a n 2 + b n 2 \frac{a_n}{\sqrt{a_n^2+b_n^2}} an2+bn2an − i b n a n 2 + b n 2 -i\frac{b_n}{\sqrt{a_n^2+b_n^2}} −ian2+bn2bn， e − i φ = c o s ( φ n ) − i ∗ s i n ( φ n ) e^{-i\varphi}=cos(\varphi_n)-i*sin(\varphi_n) e−iφ=cos(φn)−i∗sin(φn) = a n a n 2 + b n 2 \frac{a_n}{\sqrt{a_n^2+b_n^2}} an2+bn2an + i b n a n 2 + b n 2 +i\frac{b_n}{\sqrt{a_n^2+b_n^2}} +ian2+bn2bn
所以：
- a n 2 + b n 2 ∗ e i θ n \sqrt{a_n^2+b_n^2}*e^{i\theta_n} an2+bn2∗eiθn= b n + i ∗ a n b_n+i*a_n bn+i∗an
- a n 2 + b n 2 ∗ e − i θ n \sqrt{a_n^2+b_n^2}*e^{-i\theta_n} an2+bn2∗e−iθn= b n − i ∗ a n b_n-i*a_n bn−i∗an
- a n 2 + b n 2 ∗ e i φ n \sqrt{a_n^2+b_n^2}*e^{i\varphi_n} an2+bn2∗eiφn= a n − i ∗ b n a_n-i*b_n an−i∗bn
- a n 2 + b n 2 ∗ e − i φ n \sqrt{a_n^2+b_n^2}*e^{-i\varphi_n} an2+bn2∗e−iφn= a n + i ∗ b n a_n+i*b_n an+i∗bn
所以：
- a n 2 + b n 2 ∗ \sqrt{a_n^2+b_n^2}* an2+bn2∗ s i n ( n w 1 t + θ n ) sin(nw_1t+\theta_n) sin(nw1t+θn)= a n − i ∗ b n 2 e i n w 1 t + a n + i ∗ b n 2 e − i n w 1 t \frac{a_n-i*b_n}{2}e^{inw_1t}+\frac{a_n+i*b_n}{2}e^{-inw_1t} 2an−i∗bneinw1t+2an+i∗bne−inw1t
- a n 2 + b n 2 ∗ \sqrt{a_n^2+b_n^2}* an2+bn2∗ c o s ( n w 1 t + φ n ) cos(nw_1t+\varphi_n) cos(nw1t+φn)= a n − i ∗ b n 2 e i n w 1 t + a n + i ∗ b n 2 e − i n w 1 t \frac{a_n-i*b_n}{2}e^{inw_1t}+\frac{a_n+i*b_n}{2}e^{-inw_1t} 2an−i∗bneinw1t+2an+i∗bne−inw1t
所以：
- f ( t ) = a 0 + a n 2 + b n 2 ∑ n = 1 ∞ s i n ( n w 1 t + θ n ) f(t)=a_0+\sqrt{a_n^2+b_n^2}\sum\limits_{n=1}^{\infty}sin(nw_1t+\theta_n) f(t)=a0+an2+bn2n=1∑∞sin(nw1t+θn)= a 0 + a n 2 + b n 2 ∑ n = 1 ∞ c o s ( n w 1 t + φ n ) a_0+\sqrt{a_n^2+b_n^2}\sum\limits_{n=1}^{\infty}cos(nw_1t+\varphi_n) a0+an2+bn2n=1∑∞cos(nw1t+φn)= a 0 + ∑ n = 1 ∞ ( a n − i ∗ b n 2 e i n w 1 t + a n + i ∗ b n 2 e − i n w 1 t ) a_0+\sum\limits_{n=1}^{\infty}(\frac{a_n-i*b_n}{2}e^{inw_1t}+\frac{a_n+i*b_n}{2}e^{-inw_1t}) a0+n=1∑∞(2an−i∗bneinw1t+2an+i∗bne−inw1t)
又因为 b n = 2 T ∫ t 0 t 0 + T f ( t ) s i n ( n w 1 t ) d t b_n=\frac{2}{T}\int_{t_0}^{t_0+T}f(t)sin(nw_1t)d_t bn=T2∫t0t0+Tf(t)sin(nw1t)dt相对于参数n为奇函数所以 a n − i ∗ b n 2 \frac{a_n-i*b_n}{2} 2an−i∗bn( n = 1 , 2 , . . . , ∞ n=1,2,...,\infty n=1,2,...,∞)= a n + i ∗ b n 2 \frac{a_n+i*b_n}{2} 2an+i∗bn( n = − 1 , − 2 , . . . , − ∞ n=-1,-2,...,-\infty n=−1,−2,...,−∞)。现在令
- F ( n w 1 ) = a n − i ∗ b n 2 F(nw_1)=\frac{a_n-i*b_n}{2} F(nw1)=2an−i∗bn n = 1 , 2 , 3 , . . . , ∞ n=1,2,3,...,\infty n=1,2,3,...,∞， F ( 0 w 1 ) = a 0 F(0w_1)=a_0 F(0w1)=a0
则 f ( t ) = ∑ n = − ∞ ∞ F ( n w 1 ) ∗ e i n w 1 t f(t)=\sum\limits_{n=-\infty}^{\infty}F(nw_1)*e^{inw_1t} f(t)=n=−∞∑∞F(nw1)∗einw1t。现在我们也同样将时间域信号变换到了频率域。只不过这里信号 f ( t ) f(t) f(t)被分解为了一系列频率不同( w = n w 1 ， n = − ∞ , . . . , − 2 , − 1 , 0 , 1 , 2 , . . . , ∞ w=nw_1，n=-\infty,...,-2,-1,0,1,2,...,\infty w=nw1，n=−∞,...,−2,−1,0,1,2,...,∞)、幅度(复数 F ( n w 1 ) = a n − i ∗ b n 2 F(nw_1)=\frac{a_n-i*b_n}{2} F(nw1)=2an−i∗bn的模长 ， n = − ∞ , . . . , − 2 , − 1 , 0 , 1 , 2 , . . . , ∞ ，n=-\infty,...,-2,-1,0,1,2,...,\infty ，n=−∞,...,−2,−1,0,1,2,...,∞)不同和相位(复数 F ( n w 1 ) = a n − i ∗ b n 2 F(nw_1)=\frac{a_n-i*b_n}{2} F(nw1)=2an−i∗bn的相位角， n = n = − ∞ , . . . , − 2 , − 1 , 0 , 1 , 2 , . . . , ∞ n=n=-\infty,...,-2,-1,0,1,2,...,\infty n=n=−∞,...,−2,−1,0,1,2,...,∞)不同的复指数函数。这里也可以以频率( w = n w 1 ， n = 0 , 1 , 2 , . . . , w=nw_1，n=0,1,2,..., w=nw1，n=0,1,2,...,)为横坐标，幅度(复数 F ( n w 1 ) = a n − i ∗ b n 2 F(nw_1)=\frac{a_n-i*b_n}{2} F(nw1)=2an−i∗bn的模长 ， n = − ∞ , . . . , − 2 , − 1 , 0 , 1 , 2 , . . . , ∞ ，n=-\infty,...,-2,-1,0,1,2,...,\infty ，n=−∞,...,−2,−1,0,1,2,...,∞)和相位(复数 F ( n w 1 ) = a n − i ∗ b n 2 F(nw_1)=\frac{a_n-i*b_n}{2} F(nw1)=2an−i∗bn的相位角， n = − ∞ , . . . , − 2 , − 1 , 0 , 1 , 2 , . . . , ∞ n=-\infty,...,-2,-1,0,1,2,...,\infty n=−∞,...,−2,−1,0,1,2,...,∞)为纵坐标作图，这样就更显示了频率域的意义了。一个简单的示意图如图2所示（注意这里的图示不是上面的结果，只是一个简单的示意而已）。单纯从 F ( n w 1 ) F(nw_1) F(nw1)的等式中我们可以知道幅度谱的图以幅度轴偶对称，相位谱的图以幅度轴奇对称。
这里与之前不同的是出现了负频率，但是频率是个标量不存在方向一说，这里出现的负频率是因为我们引入了复指数的函数的原因，引入复指数函数完全是为了方便理论上的数学运算。这里复函数 F ( n w 1 ) = a n − i ∗ b n 2 F(nw_1)=\frac{a_n-i*b_n}{2} F(nw1)=2an−i∗bn的模长为 = a n 2 + b n 2 2 =\frac{\sqrt{a_n^2+b_n^2}}{2} =2an2+bn2，它的值为分解为正弦或余弦函数时对应的幅度值的一半(n>0)，又因为此时幅度谱的图以幅度轴偶对称，则当我们利用复指数函数分解之后得到对应频率的幅度值之后应该再乘以2才是对应的频率值的信号的真实幅度。（这一点对于实信号即 f ( x ) f(x) f(x)为实函数，我觉得应该是如此的，但是有时候似乎 f ( x ) f(x) f(x)可以为复函数即复信号，这咱也不知道了，还有复信号？）
以上的信号分解针对的是周期信号，那非周期信号又将如何处理？由以上推导可以得到
F
(
n
w
1
)
=
a
n
−
i
∗
b
n
2
=
1
T
∫
t
0
t
0
+
T
f
(
t
)
c
o
s
(
n
w
1
t
)
d
t
−
i
∗
1
T
∫
t
0
t
0
+
T
f
(
t
)
s
i
n
(
n
w
1
t
)
d
t
=
1
T
∫
t
0
t
0
+
T
f
(
t
)
(
c
o
s
(
n
w
1
t
)
−
i
∗
s
i
n
(
n
w
1
t
)
)
d
t
=
1
T
∫
t
0
t
0
+
T
f
(
t
)
e
−
i
n
w
1
t
d
t
=
1
T
∫
−
T
2
T
2
f
(
t
)
e
−
i
n
w
1
t
d
t
F(nw_1)=\frac{a_n-i*b_n}{2}=\frac{1}{T}\int_{t_0}^{t_0+T}f(t)cos(nw_1t)d_t-i*\frac{1}{T}\int_{t_0}^{t_0+T}f(t)sin(nw_1t)d_t=\frac{1}{T}\int_{t_0}^{t_0+T}f(t)(cos(nw_1t)-i*sin(nw_1t))d_t=\frac{1}{T}\int_{t_0}^{t_0+T}f(t)e^{-inw_1t}d_t=\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{-inw_1t}d_t
F(nw1)=2an−i∗bn=T1∫t0t0+Tf(t)cos(nw1t)dt−i∗T1∫t0t0+Tf(t)sin(nw1t)dt=T1∫t0t0+Tf(t)(cos(nw1t)−i∗sin(nw1t))dt=T1∫t0t0+Tf(t)e−inw1tdt=T1∫−2T2Tf(t)e−inw1tdt
n
=
1
,
2
,
.
.
.
,
∞
n=1,2,...,\infin
n=1,2,...,∞
在等式两边乘以T可以得到：
- T F ( n w 1 ) = 2 π F ( n w 1 ) w 1 = ∫ − T 2 T 2 f ( t ) e − i n w 1 t d t TF(nw_1)=\frac{2\pi F(nw_1)}{w_1}=\int_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{-inw_1t}d_t TF(nw1)=w12πF(nw1)=∫−2T2Tf(t)e−inw1tdt
当周期信号 f ( t ) f(t) f(t)的周期趋于无穷大时，频谱图的频率轴的点间间隔 w 1 w_1 w1趋向于0，离散频率 n w 1 nw_1 nw1变成了连续频率。在这种极限情况下 F ( n w 1 ) F(nw_1) F(nw1)趋向于0，但是 T F ( n w 1 ) = 2 π F ( n w 1 ) w 1 TF(nw_1)=\frac{2\pi F(nw_1)}{w_1} TF(nw1)=w12πF(nw1)可望不趋于0，而趋于有限值，且变成一个连续函数，通常记作 F ( w ) F(w) F(w)即：
- F ( w ) = lim T → ∞ 2 π F ( n w 1 ) w 1 = lim T → ∞ ∫ − T 2 T 2 f ( t ) e − i n w 1 t d t = ∫ − ∞ ∞ f ( t ) e − i w t d t F(w)=\lim\limits_{T \to \infin}\frac{2\pi F(nw_1)}{w_1}=\lim\limits_{T \to \infin}\int_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{-inw_1t}d_t=\int_{-\infin}^{\infin}f(t)e^{-iwt}d_t F(w)=T→∞limw12πF(nw1)=T→∞lim∫−2T2Tf(t)e−inw1tdt=∫−∞∞f(t)e−iwtdt（这就是非周期信号的傅里叶变换）
同理因为 f ( t ) = ∑ n = − ∞ ∞ F ( n w 1 ) ∗ e i n w 1 t = 1 2 π ∑ n = − ∞ ∞ 2 π F ( n w 1 ) w 1 ∗ e i n w 1 t ∗ w 1 f(t)=\sum\limits_{n=-\infty}^{\infty}F(nw_1)*e^{inw_1t}=\frac{1}{2\pi}\sum\limits_{n=-\infty}^{\infty}\frac{2\pi F(nw_1)}{w_1}*e^{inw_1t}*w_1 f(t)=n=−∞∑∞F(nw1)∗einw1t=2π1n=−∞∑∞w12πF(nw1)∗einw1t∗w1，当T趋向于无穷大时， w 1 w_1 w1趋向于0，由定积分的定义可得此时 f ( t ) = 1 2 π ∑ n = − ∞ ∞ F ( w ) ∗ e i w t ∗ d w = 1 2 π ∫ − ∞ ∞ F ( w ) ∗ e i w t d w f(t)=\frac{1}{2\pi}\sum\limits_{n=-\infty}^{\infty}F(w)*e^{iwt}*d_w=\frac{1}{2\pi}\int_{-\infin}^{\infin}F(w)*e^{iwt}d_w f(t)=2π1n=−∞∑∞F(w)∗eiwt∗dw=2π1∫−∞∞F(w)∗eiwtdw（这就是非周期信号的傅里叶反变换）
F ( w ) F(w) F(w)一般是复函数，可以写作 F ( w ) = ∣ F ( w ) ∣ e i φ ( w ) F(w)=|F(w)|e^{i\varphi(w)} F(w)=∣F(w)∣eiφ(w)， ∣ F ( w ) ∣ |F(w)| ∣F(w)∣是模长， φ ( w ) \varphi(w) φ(w)是相位角。与周期信号的复指数分解一样也可以画出以频率( w w w连续变量)为横坐标，幅度( ∣ F ( w ) ∣ |F(w)| ∣F(w)∣)和相位角( φ ( w ) \varphi(w) φ(w))为纵坐标的频率域图，只不过这里的图是连续的而不再是离散的点。一个简单的例子如图3所示。
可以发现这里也出现了负频率，与周期信号的复指数分解一样，出现的负频率是因为我们引入了复指数的函数的原因，引入复指数函数完全是为了方便理论上的数学运算。对于正频率w的信号的幅度值应该为计算得到的正频率w的幅度值的两倍（这一点对于实信号即 f ( x ) f(x) f(x)为实函数，我觉得应该是如此的，但是有时候似乎 f ( x ) f(x) f(x)可以为复函数即复信号，这咱也不知道了，还有复信号？）。
电子计算机是无法处理连续的信号的，一般都是对连续信号进行离散采样后再进行处理。如图4所示。
信号 f ( t ) f(t) f(t)（非周期连续信号）经过采样之后可以表示为采样信号 f s ( t ) = ∑ k = − ∞ ∞ f ( t ) ∗ δ ( t − k T ) f_s(t)=\sum\limits_{k=-\infty}^{\infty}f(t)*\delta(t-kT) fs(t)=k=−∞∑∞f(t)∗δ(t−kT)，这时得到了一系列离散的点，点与点指点的采样周期为T。但是此时 f s ( t ) f_s(t) fs(t)依然是连续信号只不过在离散点之外的其它地方值为0。 δ ( t ) \delta(t) δ(t)函数为冲击函数。对采样信号 f s ( t ) f_s(t) fs(t)作傅里叶变换可得：
- F ( w ) = ∫ − ∞ ∞ f s ( t ) e − i w t d t = ∫ − ∞ ∞ e − i w t ∑ k = − ∞ ∞ f ( t ) ∗ δ ( t − k T ) d t = ∫ − ∞ ∞ e − i w t ∑ k = − ∞ ∞ f ( k T ) ∗ δ ( t − k T ) d t = ∑ k = − ∞ ∞ f ( k T ) ∫ − ∞ ∞ e − i w t δ ( t − k T ) d t = ∑ k = − ∞ ∞ f ( k T ) e − i w k T = ∑ k = − ∞ ∞ f ( k ) e − i w k （ 令 k 表 示 k T ） F(w)=\int_{-\infin}^{\infin}f_s(t)e^{-iwt}d_t=\int_{-\infin}^{\infin}e^{-iwt}\sum\limits_{k=-\infty}^{\infty}f(t)*\delta(t-kT)d_t=\int_{-\infin}^{\infin}e^{-iwt}\sum\limits_{k=-\infty}^{\infty}f(kT)*\delta(t-kT)d_t=\sum\limits_{k=-\infty}^{\infty}f(kT)\int_{-\infin}^{\infin}e^{-iwt}\delta(t-kT)d_t=\sum\limits_{k=-\infty}^{\infty}f(kT)e^{-iwkT}=\sum\limits_{k=-\infty}^{\infty}f(k)e^{-iwk}（令k表示kT） F(w)=∫−∞∞fs(t)e−iwtdt=∫−∞∞e−iwtk=−∞∑∞f(t)∗δ(t−kT)dt=∫−∞∞e−iwtk=−∞∑∞f(kT)∗δ(t−kT)dt=k=−∞∑∞f(kT)∫−∞∞e−iwtδ(t−kT)dt=k=−∞∑∞f(kT)e−iwkT=k=−∞∑∞f(k)e−iwk（令k表示kT）（这就是非周期信号的采样信号的离散时间傅里叶变换）
因为 F ( w + 2 π ) = ∑ k = − ∞ ∞ f ( k ) e − i ( w + 2 π ) k = ∑ k = − ∞ ∞ f ( k ) e − i w k e − i 2 π k = ∑ k = − ∞ ∞ f ( k ) e − i w k = F ( w ) F(w+2\pi)=\sum\limits_{k=-\infty}^{\infty}f(k)e^{-i(w+2\pi)k}=\sum\limits_{k=-\infty}^{\infty}f(k)e^{-iwk}e^{-i2\pi k}=\sum\limits_{k=-\infty}^{\infty}f(k)e^{-iwk}=F(w) F(w+2π)=k=−∞∑∞f(k)e−i(w+2π)k=k=−∞∑∞f(k)e−iwke−i2πk=k=−∞∑∞f(k)e−iwk=F(w)因此非周期连续信号 f ( t ) f(t) f(t)的采样信号 f s ( t ) f_s(t) fs(t)的傅里叶变换是周期性的连续函数。
由傅里叶变换的反变换得 f s ( t ) = 1 2 π ∫ − ∞ ∞ F ( w ) ∗ e i w t d w = 1 2 π ∫ − ∞ ∞ ( ∑ k = − ∞ ∞ f ( k ) e − i w k ) ∗ e i w t d w = 1 2 π ∫ − ∞ ∞ ( ∑ k = − ∞ ∞ f ( k ) e i w ( t − k ) ) d w = 1 2 π ∑ k = − ∞ ∞ f ( k ) ∫ − ∞ ∞ e i w ( t − k ) ) d w = 1 2 π ∑ k = − ∞ ∞ f ( k ) ∫ − ∞ ∞ ( c o s ( w ( t − k ) ) − i s i n ( w ( t − k ) ) ) d w = 1 2 π ∑ k = − ∞ ∞ f ( k ) ∫ − ∞ ∞ ( c o s ( w ( t − k ) ) − i s i n ( w ( t − k ) ) ) d w f_s(t)=\frac{1}{2\pi}\int_{-\infin}^{\infin}F(w)*e^{iwt}d_w=\frac{1}{2\pi}\int_{-\infin}^{\infin}(\sum\limits_{k=-\infty}^{\infty}f(k)e^{-iwk})*e^{iwt}d_w=\frac{1}{2\pi}\int_{-\infin}^{\infin}(\sum\limits_{k=-\infty}^{\infty}f(k)e^{iw(t-k)})d_w=\frac{1}{2\pi}\sum\limits_{k=-\infty}^{\infty}f(k)\int_{-\infin}^{\infin}e^{iw(t-k)})d_w=\frac{1}{2\pi}\sum\limits_{k=-\infty}^{\infty}f(k)\int_{-\infin}^{\infin}(cos(w(t-k))-isin(w(t-k)))d_w=\frac{1}{2\pi}\sum\limits_{k=-\infty}^{\infty}f(k)\int_{-\infin}^{\infin}(cos(w(t-k))-isin(w(t-k)))d_w fs(t)=2π1∫−∞∞F(w)∗eiwtdw=2π1∫−∞∞(k=−∞∑∞f(k)e−iwk)∗eiwtdw=2π1∫−∞∞(k=−∞∑∞f(k)eiw(t−k))dw=2π1k=−∞∑∞f(k)∫−∞∞eiw(t−k))dw=2π1k=−∞∑∞f(k)∫−∞∞(cos(w(t−k))−isin(w(t−k)))dw=2π1k=−∞∑∞f(k)∫−∞∞(cos(w(t−k))−isin(w(t−k)))dw因为当 t ! = k t!=k t!=k时 ∫ − ∞ ∞ ( c o s ( w ( t − k ) ) − i s i n ( w ( t − k ) ) ) d w = 0 \int_{-\infin}^{\infin}(cos(w(t-k))-isin(w(t-k)))d_w=0 ∫−∞∞(cos(w(t−k))−isin(w(t−k)))dw=0所以上式化为 f s ( t ) = 1 2 π f s ( t ) ∫ − ∞ ∞ 1 d w f_s(t)=\frac{1}{2\pi}f_s(t)\int_{-\infin}^{\infin}1d_w fs(t)=2π1fs(t)∫−∞∞1dw但是到这里我们发现无穷限的积分无求了。参考资料说是这里求反变换的时候不能取无穷限因为这里F(w)是周期为 2 π 2\pi 2π的周期函数，取一个周期的积分区间足以，这时我们如果取一个周期 2 π 2\pi 2π的积分区间就可以得到反变换的公式为 f s ( t ) = 1 2 π ∫ − π π F ( w ) ∗ e i w t d w f_s(t)=\frac{1}{2\pi}\int_{-\pi}^{\pi}F(w)*e^{iwt}d_w fs(t)=2π1∫−ππF(w)∗eiwtdw（这就是非周期信号的采样信号的离散时间傅里叶反变换），这里由于我的知识储备不够我暂时也无法解释（信号与系统课本是通过Z变换导出的离散时间傅里叶变换与反变换的）
以上分析的傅里叶级数以及傅里叶变换都是针对连续信号的肯定无法用于电子计算机。对于离散时间傅里叶变换其正变换是无穷项的求和，逆变换时一个周期期间的积分这也是无法用于电子计算机取处理的，因为电子计算机肯定只能处理离散以及有限项的数据。这就引出了离散傅里叶变换。
由以上的推导可得周期信号的复指数形式的分解与合成公式为：
- F ( n w 1 ) = 1 T ∫ − T 2 T 2 f ( t ) e − i n w 1 t d t F(nw_1)=\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{-inw_1t}d_t F(nw1)=T1∫−2T2Tf(t)e−inw1tdt n = , − ∞ , . . . , − 2 , − 1 , 0 , 1 , 2 , . . . , ∞ n=,-\infin,...,-2,-1,0,1,2,...,\infin n=,−∞,...,−2,−1,0,1,2,...,∞
- f ( t ) = ∑ n = − ∞ ∞ F ( n w 1 ) ∗ e i n w 1 t f(t)=\sum\limits_{n=-\infty}^{\infty}F(nw_1)*e^{inw_1t} f(t)=n=−∞∑∞F(nw1)∗einw1t
和离散时间傅里叶变换一样，周期连续信号 f ( t ) f(t) f(t)（周期为T）经过采样之后可以表示为采样信号 f s ( t ) = ∑ k = − ∞ ∞ f ( t ) ∗ δ ( t − k Δ t ) f_s(t)=\sum\limits_{k=-\infty}^{\infty}f(t)*\delta(t-k\Delta{t}) fs(t)=k=−∞∑∞f(t)∗δ(t−kΔt)，这时得到了一系列离散的点，点与点之间的采样周期为 Δ t \Delta{t} Δt，采样频率 f = 1 Δ t f=\frac{1}{\Delta{t}} f=Δt1。但是此时 f s ( t ) f_s(t) fs(t)依然是连续周期信号只不过在离散点之外的其它地方值为0。这时一个周期采样N个点，如图5所示。则采样信号 f s ( t ) f_s(t) fs(t)的复指数傅里叶级数的幅度函数为： F ( n w 1 ) = 1 T ∫ 0 T f ( t ) e − i n w 1 t d t F(nw_1)=\frac{1}{T}\int_{0}^{T}f(t)e^{-inw_1t}d_t F(nw1)=T1∫0Tf(t)e−inw1tdt。因为采样后采样信号只在N个点处有值其它地方均为0， N Δ t = T N\Delta{t}=T NΔt=T则上式可以化简为 F ( n w 1 ) = Δ t T ∑ k = 0 N − 1 f ( k Δ t ) e − i n w 1 k Δ t = 1 N ∑ k = 0 N − 1 f ( k Δ t ) e − i 2 π k n N F(nw_1)=\frac{\Delta{t}}{T}\sum\limits_{k=0}^{N-1}f(k\Delta{t})e^{-inw_1k\Delta{t}}=\frac{1}{N}\sum\limits_{k=0}^{N-1}f(k\Delta{t})e^{\frac{-i2\pi kn}{N}} F(nw1)=TΔtk=0∑N−1f(kΔt)e−inw1kΔt=N1k=0∑N−1f(kΔt)eN−i2πkn再令 f ( k Δ t ) = f ( k ) f(k\Delta{t})=f(k) f(kΔt)=f(k)， F ( n w 1 ) = F ( n ) F(nw_1)=F(n) F(nw1)=F(n)，则：
- F ( n ) = 1 N ∑ k = 0 N − 1 f ( k ) e − i 2 π k n N F(n)=\frac{1}{N}\sum\limits_{k=0}^{N-1}f(k)e^{\frac{-i2\pi kn}{N}} F(n)=N1k=0∑N−1f(k)eN−i2πkn（正变换）
这时我们可以发现 F ( n ) F(n) F(n)是以N为周期的周期函数。和离散时间傅里叶变换的推导一样，由于这里 F ( n ) F(n) F(n)为周期函数，我们将 F ( n ) F(n) F(n)带入周期信号的复指数分解的公式中时求和范围应该是一个周期，同时采样信号 f s ( t ) f_s(t) fs(t)只在采样点处有值，其它地方的值为0，因此：
- f s ( t ) = f s ( m Δ t ) = f ( m ) = ∑ n = 0 N − 1 e i n w 1 m Δ t 1 N ∑ k = 0 N − 1 f ( k ) e − i 2 π k n N = 1 N ∑ n = 0 N − 1 e i n 2 π N Δ t m Δ t ∑ k = 0 N − 1 f ( k ) e − i 2 π k n N = 1 N ∑ n = 0 N − 1 e i n 2 π N m ∑ k = 0 N − 1 f ( k ) e − i 2 π k n N f_s(t)=f_s(m\Delta t)=f(m)=\sum\limits_{n=0}^{N-1}e^{inw_1m\Delta t}\frac{1}{N}\sum\limits_{k=0}^{N-1}f(k)e^{\frac{-i2\pi kn}{N}}=\frac{1}{N}\sum\limits_{n=0}^{N-1}e^{in\frac{2\pi}{N\Delta t}m\Delta t}\sum\limits_{k=0}^{N-1}f(k)e^{\frac{-i2\pi kn}{N}}=\frac{1}{N}\sum\limits_{n=0}^{N-1}e^{in\frac{2\pi}{N}m}\sum\limits_{k=0}^{N-1}f(k)e^{\frac{-i2\pi kn}{N}} fs(t)=fs(mΔt)=f(m)=n=0∑N−1einw1mΔtN1k=0∑N−1f(k)eN−i2πkn=N1n=0∑N−1einNΔt2πmΔtk=0∑N−1f(k)eN−i2πkn=N1n=0∑N−1einN2πmk=0∑N−1f(k)eN−i2πkn
现在考虑：
- S N = ∑ k = 0 N − 1 e i k α = ∑ k = 0 N − 1 ( e i α ) k = 1 − ( e i α ) N 1 − e i α S_N=\sum\limits_{k=0}^{N-1}e^{ik\alpha}=\sum\limits_{k=0}^{N-1}{({e^{i\alpha})}^k}=\frac{1-({e^{i\alpha})}^N}{1-e^{i\alpha}} SN=k=0∑N−1eikα=k=0∑N−1(eiα)k=1−eiα1−(eiα)N，当 α = 0 , ± 2 π , ± 4 π , . . . \alpha=0,\pm2\pi,\pm4\pi,... α=0,±2π,±4π,...时， S N = N S_N=N SN=N
如果希望当 α ! = 0 , ± 2 π , ± 4 π , . . . \alpha!=0,\pm2\pi,\pm4\pi,... α!=0,±2π,±4π,...时 S N = 0 S_N=0 SN=0，则 e i α ！ = 1 e^{i\alpha}！=1 eiα！=1， e i α N = 1 e^{i\alpha N}=1 eiαN=1，即 α = 2 π n N \alpha=\frac{2\pi n}{N} α=N2πn，n为整数且 n ！ = 0 ， ± N , ± 2 N , . . . . n！=0，\pm N,\pm 2N,.... n！=0，±N,±2N,....，因此：
S N = ∑ k = 0 N − 1 e i k 2 π n N = { N n = 0 ， ± N , ± 2 N , . . . . 0 n ！ = 0 ， ± N , ± 2 N , . . . . S_N=\sum\limits_{k=0}^{N-1}e^{ik\frac{2\pi n}{N}} =\begin{cases} N& \text{$n=0，\pm N,\pm 2N,....$}\\ 0& \text{$n！=0，\pm N,\pm 2N,....$} \end{cases} SN=k=0∑N−1eikN2πn={N0n=0，±N,±2N,....n！=0，±N,±2N,....
因此 f s ( t ) = f s ( m Δ t ) = f ( m ) = 1 N ∑ n = 0 N − 1 e i n 2 π N m ∑ k = 0 N − 1 f ( k ) e − i 2 π k n N = 1 N ∑ n = 0 N − 1 ∑ k = 0 N − 1 f ( k ) e i 2 π ( m − k ) n N = 1 N ∑ k = 0 N − 1 f ( k ) ∑ n = 0 N − 1 e i 2 π ( m − k ) n N f_s(t)=f_s(m\Delta t)=f(m)=\frac{1}{N}\sum\limits_{n=0}^{N-1}e^{in\frac{2\pi}{N}m}\sum\limits_{k=0}^{N-1}f(k)e^{\frac{-i2\pi kn}{N}}=\frac{1}{N}\sum\limits_{n=0}^{N-1}\sum\limits_{k=0}^{N-1}f(k)e^{\frac{i2\pi (m-k)n}{N}}=\frac{1}{N}\sum\limits_{k=0}^{N-1}f(k)\sum\limits_{n=0}^{N-1}e^{\frac{i2\pi (m-k)n}{N}} fs(t)=fs(mΔt)=f(m)=N1n=0∑N−1einN2πmk=0∑N−1f(k)eN−i2πkn=N1n=0∑N−1k=0∑N−1f(k)eNi2π(m−k)n=N1k=0∑N−1f(k)n=0∑N−1eNi2π(m−k)n，所以当 m − k = 0 ， ± N , ± 2 N , . . . . m-k=0，\pm N,\pm 2N,.... m−k=0，±N,±2N,....时 f s ( t ) = f s ( m Δ t ) = f ( m ) = 1 N ∗ N ∗ f ( k ) = f ( k ) f_s(t)=f_s(m\Delta t)=f(m)=\frac{1}{N}*N*f(k)=f(k) fs(t)=fs(mΔt)=f(m)=N1∗N∗f(k)=f(k)（逆变换）。因为这里 f ( t ) f(t) f(t)和 F ( n ) F(n) F(n)都是周期函数可将 f s ( t ) f_s(t) fs(t)和 F ( n ) F(n) F(n)都限制在一个周期内，因此将m和n的范围限制在 0 到 N − 1 0到N-1 0到N−1。此时正变换和逆变换都将是求和公式而不再是积分且求和也只在有限个值N。这样就可以应用在电子计算中。现在给出离散傅里叶变换的正反变换公式：
- F ( n ) = ∑ k = 0 N − 1 f ( k ) e − i 2 π k n N F(n)=\sum\limits_{k=0}^{N-1}f(k)e^{\frac{-i2\pi kn}{N}} F(n)=k=0∑N−1f(k)eN−i2πkn（离散傅里叶变换正变换）
- f ( m ) = 1 N ∑ n = 0 N − 1 F ( n ) e i n 2 π N m f(m)=\frac{1}{N}\sum\limits_{n=0}^{N-1}F(n)e^{in\frac{2\pi}{N}m} f(m)=N1n=0∑N−1F(n)einN2πm（离散傅里叶变换逆变换）
注意上面将 1 N \frac{1}{N} N1挪到了逆变换公式的前面，但是有的资料中没有移动还是在正变换公式中。同时 F ( N − n ) = ∑ k = 0 N − 1 f ( k ) e − i 2 π k ( N − n ) N = ∑ k = 0 N − 1 f ( k ) e i 2 π k n N F(N-n)=\sum\limits_{k=0}^{N-1}f(k)e^{\frac{-i2\pi k(N-n)}{N}}=\sum\limits_{k=0}^{N-1}f(k)e^{\frac{i2\pi kn}{N}} F(N−n)=k=0∑N−1f(k)eN−i2πk(N−n)=k=0∑N−1f(k)eNi2πkn，因此可知 F ( n ) F(n) F(n)和 F ( N − n ) F(N-n) F(N−n)实部相同虚部相反，但是 F ( n ) F(n) F(n)和 F ( N − n ) F(N-n) F(N−n)的模长即幅度还是相同的。当 N = 8 N=8 N=8时（在离散傅里叶变换中为了方便计算机处理N的值一般取为2的幂）， F ( n ) F(n) F(n)的一个简单的实例如图6所示。
图6中 n = 0 , ± 8 n=0,\pm 8 n=0,±8处的点是一个周期的开始也是一个周期的结束。因为采样频率 f = 1 Δ t f=\frac{1}{\Delta t} f=Δt1，周期信号 f ( t ) f(t) f(t)的周期 T = 8 Δ t T=8\Delta t T=8Δt，点 n , n = 0 , 1 , 2 , 3 , 4 , 5 , 6 , 7 n,n=0,1,2,3,4,5,6,7 n,n=0,1,2,3,4,5,6,7处对应的频率值为 n 1 8 Δ t = n 8 ∗ f n\frac{1}{8\Delta t}=\frac{n}{8}*f n8Δt1=8n∗f。由采样定理可知要从采样信号中恢复原信号则采样频率必须大于等于该信号中最高频率分量的两倍。因此如果原信号中有大于 f 2 \frac{f}{2} 2f的频率分量该频率分量就无法恢复了（对于离散傅里叶变换就会出现频谱泄露还有啥频谱混叠之类的问题，由于知识储备不够这里就不解释了），因此一般在得到傅里叶变换的结果之后我们只取对应于点0到点 N 2 \frac{N}{2} 2N之间的数据。从图6中可以看到也有负频率存在，就如之前提到的这时由于引入复指数函数的原因。因此对应频率点的实际幅度值应该为计算得到的两倍，但是不知道为啥我看到的一些代码中最后没有进行乘2的操作。
最后来讲一下快速傅里叶变换，快速傅里叶变换其实没啥。对于给定的N个采样点数据，离散傅里叶变换和快速傅里叶变换算出的结果一样。但是快速傅里叶变换比离散傅里叶变换要快很对特别是当N的值变大时。也正是由于快，使得一些应用称为现实。如果一个应用需要在10秒内完成傅里叶变换，结果离散傅里叶变换要算一个小时，那就没啥用了。下面具体讲解傅里叶变换。
在以下推导过程中令 W N = e − j 2 π N W_N=e^{{-j}{\frac{2\pi} {N}}} WN=e−jN2π， { W N } k = W N k \{W_N\}^k=W_N^k {WN}k=WNk且时间域信号采样数据点 f [ k ] f[k] f[k]和频率域数据点 F [ n ] F[n] F[n]均为复数，则可以得到（j为复数单位）：
- W N 0 = − W N N 2 W_N^0=-W_N^{\frac{N} {2}} WN0=−WN2N， W N 1 = − W N N 2 + 1 W_N^1=-W_N^{\frac{N} {2}+1} WN1=−WN2N+1，…， W N N 2 − 1 = − W N N 2 + N 2 − 1 = − W N N − 1 W_N^{\frac{N} {2}-1}=-W_N^{\frac{N}{2}+\frac{N}{2}-1}=-W_N^{N-1} WN2N−1=−WN2N+2N−1=−WNN−1，也即 W N i = − W N N 2 + i W_N^{i}=-W_N^{\frac{N}{2}+i} WNi=−WN2N+i。（公式一）
- W N m ∗ ( n + N ) = W N m ∗ n W_N^{m*(n+N)}=W_N^{m*n} WNm∗(n+N)=WNm∗n。（公式二）
因为
W
N
N
2
+
i
W_N^{\frac{N}{2}+i}
WN2N+i=
e
−
j
2
π
N
(
N
2
+
i
)
e^{{-j}{\frac{2\pi} {N}}({\frac{N}{2}+i})}
e−jN2π(2N+i)=
e
−
j
π
∗
e
−
j
2
π
N
i
e^{{-j}\pi}*e^{{-j}{\frac{2\pi} {N}}i}
e−jπ∗e−jN2πi=
W
N
i
∗
{
c
o
s
(
π
)
−
j
∗
s
i
n
(
π
)
}
W_N^{i}*\{cos(\pi)-j*sin(\pi)\}
WNi∗{cos(π)−j∗sin(π)}=
−
W
N
i
-W_N^{i}
−WNi因此公式一得证。
W
N
m
∗
(
n
+
N
)
W_N^{m*(n+N)}
WNm∗(n+N)=
e
−
j
2
π
N
m
(
n
+
N
)
e^{{-j}{\frac{2\pi} {N}}m(n+N)}
e−jN2πm(n+N)=
e
−
j
2
π
N
m
n
e^{{-j}{\frac{2\pi} {N}}mn}
e−jN2πmn
∗
*
∗
e
−
j
2
π
N
m
N
e^{{-j}{\frac{2\pi} {N}}mN}
e−jN2πmN=
W
N
m
∗
n
W_N^{m*n}
WNm∗n
∗
*
∗
e
−
j
2
π
N
m
N
e^{{-j}{\frac{2\pi} {N}}mN}
e−jN2πmN=
W
N
m
∗
n
W_N^{m*n}
WNm∗n
∗
*
∗
e
−
j
2
π
m
e^{{-j}{{2\pi}}m}
e−j2πm=
W
N
m
∗
n
∗
{
c
o
s
(
m
∗
2
π
)
−
j
∗
s
i
n
(
m
∗
2
π
)
}
W_N^{m*n}*\{cos(m*2\pi)-j*sin(m*2\pi)\}
WNm∗n∗{cos(m∗2π)−j∗sin(m∗2π)}=
W
N
m
∗
n
W_N^{m*n}
WNm∗n因此公式二得证。
离散傅里叶变换的公式如下：
- F [ n ] = ∑ k = 0 N − 1 f [ k ] e − j 2 π N k n F[n]=\sum\limits_{k=0}^{N-1}f[k]e^{{-j}{\frac{2\pi}{N}}kn} F[n]=k=0∑N−1f[k]e−jN2πkn= ∑ k = 0 N − 1 f [ k ] W N k n \sum\limits_{k=0}^{N-1}f[k]W_N^{kn} k=0∑N−1f[k]WNkn= ∑ m = 0 N 2 − 1 f [ 2 m ] W N 2 m n + ∑ m = 0 N 2 − 1 f [ 2 m + 1 ] W N ( 2 m + 1 ) n \sum\limits_{m=0}^{{\frac{N}{2}}-1}f[2m]W_N^{2mn}+\sum\limits_{m=0}^{{\frac{N}{2}}-1}f[2m+1]W_N^{(2m+1)n} m=0∑2N−1f[2m]WN2mn+m=0∑2N−1f[2m+1]WN(2m+1)n= ∑ m = 0 N 2 − 1 f [ 2 m ] W N 2 m n + W N n ∑ m = 0 N 2 − 1 f [ 2 m + 1 ] W N 2 m n \sum\limits_{m=0}^{{\frac{N}{2}}-1}f[2m]W_{{\frac{N}{2}}}^{mn}+W_N^{n}\sum\limits_{m=0}^{{\frac{N}{2}}-1}f[2m+1]W_{{\frac{N}{2}}}^{mn} m=0∑2N−1f[2m]W2Nmn+WNnm=0∑2N−1f[2m+1]W2Nmn
由此可知在N个时间域采样点上计算频率域信号的数据点时可以将时间域采样点分作两部分分别计算，一部分在索引为偶数的时间域采样点上先进行计算，一部分在索引为奇数的时间域采样点上先进行计算，最后通过这两部分来得到最后的结果。当时间域的信号的采样点的个数为8时（一般在进行快速傅里叶变换时都会要求进行计算的时间域采样点的个数为2的幂，这样做是为了方便计算），以上公式可以化为：
- F [ n ] = ∑ k = 0 7 f [ k ] e − j 2 π 8 k n F[n]=\sum\limits_{k=0}^{7}f[k]e^{{-j}{\frac{2\pi} {8}}kn} F[n]=k=0∑7f[k]e−j82πkn= ∑ m = 0 3 f [ 2 m ] W 4 m n + W 8 n ∑ m = 0 3 f [ 2 m + 1 ] W 4 m n \sum\limits_{m=0}^{3}f[2m]W_{4}^{mn}+W_8^{n}\sum\limits_{m=0}^{3}f[2m+1]W_{4}^{mn} m=0∑3f[2m]W4mn+W8nm=0∑3f[2m+1]W4mn
现在令 b [ n ] = ∑ m = 0 3 f [ 2 m ] W 4 m n b[n]=\sum\limits_{m=0}^{3}f[2m]W_{4}^{mn} b[n]=m=0∑3f[2m]W4mn， c [ n ] = ∑ m = 0 3 f [ 2 m + 1 ] W 4 m n c[n]=\sum\limits_{m=0}^{3}f[2m+1]W_{4}^{mn} c[n]=m=0∑3f[2m+1]W4mn， n = 0 , 1 , 2 , 3 , 4 , 5 , 6 , 7 n=0,1,2,3,4,5,6,7 n=0,1,2,3,4,5,6,7，因此可以得到：
- F [ 0 ] = b [ 0 ] + W 8 0 ∗ c [ 0 ] F[0]=b[0]+W_8^{0}*c[0] F[0]=b[0]+W80∗c[0]
- F [ 1 ] = b [ 1 ] + W 8 1 ∗ c [ 1 ] F[1]=b[1]+W_8^{1}*c[1] F[1]=b[1]+W81∗c[1]
- F [ 2 ] = b [ 2 ] + W 8 2 ∗ c [ 2 ] F[2]=b[2]+W_8^{2}*c[2] F[2]=b[2]+W82∗c[2]
- F [ 3 ] = b [ 3 ] + W 8 3 ∗ c [ 3 ] F[3]=b[3]+W_8^{3}*c[3] F[3]=b[3]+W83∗c[3]
- F [ 4 ] = b [ 4 ] + W 8 4 ∗ c [ 4 ] F[4]=b[4]+W_8^{4}*c[4] F[4]=b[4]+W84∗c[4]
- F [ 5 ] = b [ 5 ] + W 8 5 ∗ c [ 5 ] F[5]=b[5]+W_8^{5}*c[5] F[5]=b[5]+W85∗c[5]
- F [ 6 ] = b [ 6 ] + W 8 6 ∗ c [ 6 ] F[6]=b[6]+W_8^{6}*c[6] F[6]=b[6]+W86∗c[6]
- F [ 7 ] = b [ 7 ] + W 8 7 ∗ c [ 7 ] F[7]=b[7]+W_8^{7}*c[7] F[7]=b[7]+W87∗c[7]
通过公式一和公式二，上面的等式可以化为：
- F [ 0 ] = b [ 0 ] + W 8 0 ∗ c [ 0 ] F[0]=b[0]+W_8^{0}*c[0] F[0]=b[0]+W80∗c[0]
- F [ 1 ] = b [ 1 ] + W 8 1 ∗ c [ 1 ] F[1]=b[1]+W_8^{1}*c[1] F[1]=b[1]+W81∗c[1]
- F [ 2 ] = b [ 2 ] + W 8 2 ∗ c [ 2 ] F[2]=b[2]+W_8^{2}*c[2] F[2]=b[2]+W82∗c[2]
- F [ 3 ] = b [ 3 ] + W 8 3 ∗ c [ 3 ] F[3]=b[3]+W_8^{3}*c[3] F[3]=b[3]+W83∗c[3]
- F [ 0 ] = b [ 0 ] − W 8 0 ∗ c [ 0 ] F[0]=b[0]-W_8^{0}*c[0] F[0]=b[0]−W80∗c[0]
- F [ 1 ] = b [ 1 ] − W 8 1 ∗ c [ 1 ] F[1]=b[1]-W_8^{1}*c[1] F[1]=b[1]−W81∗c[1]
- F [ 2 ] = b [ 2 ] − W 8 2 ∗ c [ 2 ] F[2]=b[2]-W_8^{2}*c[2] F[2]=b[2]−W82∗c[2]
- F [ 3 ] = b [ 3 ] − W 8 3 ∗ c [ 3 ] F[3]=b[3]-W_8^{3}*c[3] F[3]=b[3]−W83∗c[3]
通过简化以后的等式可知，这样操作可以减少一半的计算量，在得到前 N 2 \frac{N}2 2N的频率域的点的值之后，只需要变换公式中的一个符号就可以得到后 N 2 \frac{N}2 2N的频率域点的值。以上等式用图形表示如图1所示。
同理可以得到：
- b [ n ] = ∑ m = 0 3 f [ 2 m ] W 4 m n b[n]=\sum\limits_{m=0}^{3}f[2m]W_{4}^{mn} b[n]=m=0∑3f[2m]W4mn= ∑ l = 0 1 f [ 4 l ] W 2 l n + W 4 n ∑ l = 0 1 f [ 4 l + 2 ] W 2 l n \sum\limits_{l=0}^{1}f[4l]W_{2}^{ln}+W_4^{n}\sum\limits_{l=0}^{1}f[4l+2]W_{2}^{ln} l=0∑1f[4l]W2ln+W4nl=0∑1f[4l+2]W2ln， n = 0 , 1 , 2 , 3 n=0,1,2,3 n=0,1,2,3
- c [ n ] = ∑ m = 0 3 f [ 2 m + 1 ] W 4 m n c[n]=\sum\limits_{m=0}^{3}f[2m+1]W_{4}^{mn} c[n]=m=0∑3f[2m+1]W4mn= ∑ s = 0 1 f [ 4 s + 1 ] W 2 s n + W 4 n ∑ s = 0 1 f [ 4 s + 3 ] W 2 s n \sum\limits_{s=0}^{1}f[4s+1]W_{2}^{sn}+W_4^{n}\sum\limits_{s=0}^{1}f[4s+3]W_{2}^{sn} s=0∑1f[4s+1]W2sn+W4ns=0∑1f[4s+3]W2sn， n = 0 , 1 , 2 , 3 n=0,1,2,3 n=0,1,2,3
现在令:
- d [ n ] = ∑ l = 0 1 f [ 4 l ] W 2 l n d[n]=\sum\limits_{l=0}^{1}f[4l]W_{2}^{ln} d[n]=l=0∑1f[4l]W2ln， e [ n ] = ∑ l = 0 1 f [ 4 l + 2 ] W 2 l n e[n]=\sum\limits_{l=0}^{1}f[4l+2]W_{2}^{ln} e[n]=l=0∑1f[4l+2]W2ln， n = 0 , 1 , 2 , 3 n=0,1,2,3 n=0,1,2,3
- f [ n ] = ∑ s = 0 1 f [ 4 s + 1 ] W 2 s n f[n]=\sum\limits_{s=0}^{1}f[4s+1]W_{2}^{sn} f[n]=s=0∑1f[4s+1]W2sn， g [ n ] = ∑ s = 0 1 f [ 4 s + 3 ] W 2 s n g[n]=\sum\limits_{s=0}^{1}f[4s+3]W_{2}^{sn} g[n]=s=0∑1f[4s+3]W2sn， n = 0 , 1 , 2 , 3 n=0,1,2,3 n=0,1,2,3
则可以得到：
- b [ 0 ] = d [ 0 ] + W 4 0 ∗ e [ 0 ] b[0]=d[0]+W_4^{0}*e[0] b[0]=d[0]+W40∗e[0]
- b [ 1 ] = d [ 1 ] + W 4 1 ∗ e [ 1 ] b[1]=d[1]+W_4^{1}*e[1] b[1]=d[1]+W41∗e[1]
- b [ 2 ] = d [ 2 ] + W 4 2 ∗ e [ 2 ] b[2]=d[2]+W_4^{2}*e[2] b[2]=d[2]+W42∗e[2]
-
b
[
3
]
=
d
[
3
]
+
W
4
3
∗
e
[
3
]
b[3]=d[3]+W_4^{3}*e[3]
b[3]=d[3]+W43∗e[3]
、、、、、、、、、、、、、、、、、、 - c [ 0 ] = f [ 0 ] + W 4 0 ∗ g [ 0 ] c[0]=f[0]+W_4^{0}*g[0] c[0]=f[0]+W40∗g[0]
- c [ 1 ] = f [ 1 ] + W 4 1 ∗ g [ 1 ] c[1]=f[1]+W_4^{1}*g[1] c[1]=f[1]+W41∗g[1]
- c [ 2 ] = f [ 2 ] + W 4 2 ∗ g [ 2 ] c[2]=f[2]+W_4^{2}*g[2] c[2]=f[2]+W42∗g[2]
- c [ 3 ] = f [ 3 ] + W 4 3 ∗ g [ 3 ] c[3]=f[3]+W_4^{3}*g[3] c[3]=f[3]+W43∗g[3]
通过公式一和公式二，上面的等式可以化为：
- b [ 0 ] = d [ 0 ] + W 4 0 ∗ e [ 0 ] b[0]=d[0]+W_4^{0}*e[0] b[0]=d[0]+W40∗e[0]
- b [ 1 ] = d [ 1 ] + W 4 1 ∗ e [ 1 ] b[1]=d[1]+W_4^{1}*e[1] b[1]=d[1]+W41∗e[1]
- b [ 2 ] = d [ 0 ] − W 4 0 ∗ e [ 0 ] b[2]=d[0]-W_4^{0}*e[0] b[2]=d[0]−W40∗e[0]
-
b
[
3
]
=
d
[
1
]
−
W
4
1
∗
e
[
1
]
b[3]=d[1]-W_4^{1}*e[1]
b[3]=d[1]−W41∗e[1]
、、、、、、、、、、、、、、、、、 - c [ 0 ] = f [ 0 ] + W 4 0 ∗ g [ 0 ] c[0]=f[0]+W_4^{0}*g[0] c[0]=f[0]+W40∗g[0]
- c [ 1 ] = f [ 1 ] + W 4 1 ∗ g [ 1 ] c[1]=f[1]+W_4^{1}*g[1] c[1]=f[1]+W41∗g[1]
- c [ 2 ] = f [ 0 ] − W 4 0 ∗ g [ 0 ] c[2]=f[0]-W_4^{0}*g[0] c[2]=f[0]−W40∗g[0]
- c [ 3 ] = f [ 1 ] − W 4 1 ∗ g [ 1 ] c[3]=f[1]-W_4^{1}*g[1] c[3]=f[1]−W41∗g[1]
以上等式用图形表示如图2所示。
同理可以得到：
- d [ n ] = ∑ l = 0 1 f [ 4 l ] W 2 l n d[n]=\sum\limits_{l=0}^{1}f[4l]W_{2}^{ln} d[n]=l=0∑1f[4l]W2ln= ∑ u = 0 0 f [ 8 u ] W 1 u n + W 2 n ∑ u = 0 0 f [ 8 u + 4 ] W 1 u n \sum\limits_{u=0}^{0}f[8u]W_{1}^{un}+W_2^{n}\sum\limits_{u=0}^{0}f[8u+4]W_{1}^{un} u=0∑0f[8u]W1un+W2nu=0∑0f[8u+4]W1un= f [ 0 ] + W 2 n f [ 4 ] f[0]+W_{2}^{n}f[4] f[0]+W2nf[4]， n = 0 , 1 n=0,1 n=0,1
- e [ n ] = ∑ l = 0 1 f [ 4 l + 2 ] W 2 l n e[n]=\sum\limits_{l=0}^{1}f[4l+2]W_{2}^{ln} e[n]=l=0∑1f[4l+2]W2ln= ∑ v = 0 0 f [ 8 v + 2 ] W 1 v n + W 2 n ∑ v = 0 0 f [ 8 v + 6 ] W 1 v n \sum\limits_{v=0}^{0}f[8v+2]W_{1}^{vn}+W_2^{n}\sum\limits_{v=0}^{0}f[8v+6]W_{1}^{vn} v=0∑0f[8v+2]W1vn+W2nv=0∑0f[8v+6]W1vn= f [ 2 ] + W 2 n f [ 6 ] f[2]+W_{2}^{n}f[6] f[2]+W2nf[6]， n = 0 , 1 n=0,1 n=0,1
- f [ n ] = ∑ s = 0 1 f [ 4 s + 1 ] W 2 s n f[n]=\sum\limits_{s=0}^{1}f[4s+1]W_{2}^{sn} f[n]=s=0∑1f[4s+1]W2sn= ∑ w = 0 0 f [ 8 w + 1 ] W 1 w n + W 2 n ∑ w = 0 1 f [ 8 w + 5 ] W 1 w n \sum\limits_{w=0}^{0}f[8w+1]W_{1}^{wn}+W_2^{n}\sum\limits_{w=0}^{1}f[8w+5]W_{1}^{wn} w=0∑0f[8w+1]W1wn+W2nw=0∑1f[8w+5]W1wn= f [ 1 ] + W 2 n f [ 5 ] f[1]+W_{2}^{n}f[5] f[1]+W2nf[5]， n = 0 , 1 n=0,1 n=0,1
- g [ n ] = ∑ s = 0 1 f [ 4 s + 3 ] W 2 s n g[n]=\sum\limits_{s=0}^{1}f[4s+3]W_{2}^{sn} g[n]=s=0∑1f[4s+3]W2sn= ∑ x = 0 0 f [ 8 x + 3 ] W 1 x n + W 2 n ∑ x = 0 0 f [ 8 x + 7 ] W 1 x n \sum\limits_{x=0}^{0}f[8x+3]W_{1}^{xn}+W_2^{n}\sum\limits_{x=0}^{0}f[8x+7]W_{1}^{xn} x=0∑0f[8x+3]W1xn+W2nx=0∑0f[8x+7]W1xn= f [ 3 ] + W 2 n f [ 7 ] f[3]+W_{2}^{n}f[7] f[3]+W2nf[7]， n = 0 , 1 n=0,1 n=0,1
化简后可得：
- d [ 0 ] = f [ 0 ] + W 2 0 f [ 4 ] d[0]=f[0]+W_{2}^{0}f[4] d[0]=f[0]+W20f[4]= f [ 0 ] + W 2 0 f [ 4 ] f[0]+W_{2}^{0}f[4] f[0]+W20f[4]= f [ 0 ] + 1 ∗ f [ 4 ] f[0]+1*f[4] f[0]+1∗f[4]
-
d
[
1
]
=
f
[
0
]
+
W
2
1
f
[
4
]
d[1]=f[0]+W_{2}^{1}f[4]
d[1]=f[0]+W21f[4]=
f
[
0
]
−
W
2
0
f
[
4
]
f[0]-W_{2}^{0}f[4]
f[0]−W20f[4]=
f
[
0
]
−
1
∗
f
[
4
]
f[0]-1*f[4]
f[0]−1∗f[4]
、、、、、、、、、、、、、、、、、、、、、、、、、、、、、 - e [ 0 ] = f [ 2 ] + W 2 0 f [ 6 ] e[0]=f[2]+W_{2}^{0}f[6] e[0]=f[2]+W20f[6]= f [ 2 ] + W 2 0 f [ 6 ] f[2]+W_{2}^{0}f[6] f[2]+W20f[6]= f [ 2 ] + 1 ∗ f [ 6 ] f[2]+1*f[6] f[2]+1∗f[6]
-
e
[
1
]
=
f
[
2
]
+
W
2
1
f
[
6
]
e[1]=f[2]+W_{2}^{1}f[6]
e[1]=f[2]+W21f[6]=
f
[
2
]
−
W
2
0
f
[
6
]
f[2]-W_{2}^{0}f[6]
f[2]−W20f[6]=
f
[
2
]
−
1
∗
f
[
6
]
f[2]-1*f[6]
f[2]−1∗f[6]
、、、、、、、、、、、、、、、、、、、、、、、、、、、、、 - f [ 0 ] = f [ 1 ] + W 2 0 f [ 5 ] f[0]=f[1]+W_{2}^{0}f[5] f[0]=f[1]+W20f[5]= f [ 1 ] + W 2 0 f [ 5 ] f[1]+W_{2}^{0}f[5] f[1]+W20f[5]= f [ 1 ] + 1 ∗ f [ 5 ] f[1]+1*f[5] f[1]+1∗f[5]
-
f
[
1
]
=
f
[
1
]
+
W
2
1
f
[
5
]
f[1]=f[1]+W_{2}^{1}f[5]
f[1]=f[1]+W21f[5]=
f
[
1
]
−
W
2
0
f
[
5
]
f[1]-W_{2}^{0}f[5]
f[1]−W20f[5]=
f
[
1
]
−
1
∗
f
[
5
]
f[1]-1*f[5]
f[1]−1∗f[5]
、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、 - g [ 0 ] = f [ 3 ] + W 2 0 f [ 7 ] g[0]=f[3]+W_{2}^{0}f[7] g[0]=f[3]+W20f[7]= f [ 3 ] + W 2 0 f [ 7 ] f[3]+W_{2}^{0}f[7] f[3]+W20f[7]= f [ 3 ] + 1 ∗ f [ 7 ] f[3]+1*f[7] f[3]+1∗f[7]
- g [ 1 ] = f [ 3 ] + W 2 1 f [ 7 ] g[1]=f[3]+W_{2}^{1}f[7] g[1]=f[3]+W21f[7]= f [ 3 ] − W 2 0 f [ 7 ] f[3]-W_{2}^{0}f[7] f[3]−W20f[7]= f [ 3 ] − 1 ∗ f [ 7 ] f[3]-1*f[7] f[3]−1∗f[7]
以上等式用图形表示如图3所示。
以上 log 2 8 = 3 \log_28=3 log28=3步的整个流程图如图4所示。
经过以上例子可以归纳出快速傅里叶算法的一般步骤为：
- 对N个时间域信号采样点数据重新排序(信号采样点索引值从0到 N − 1 N-1 N−1)，使得采样点 f [ k ] f[k] f[k]的数据与采样点 f [ k B i t R e v e r s e ] f[kBitReverse] f[kBitReverse]的数据交换。kBitReverse的值是K的二进制值的比特位逆序之后的值。对于有 N = 8 N=8 N=8个时间域信号采样点的重排序例子图5所示。
- 对N个时间域信号采样点数据重新排好序后需要进过 log 2 N \log_2N log2N步的迭代才能得到最后的频率域信号，每一步的迭代都会计算子DFT，每个子DFT的计算由若干个蝴蝶(butterfly)计算单元组成,如图6所示。
快速傅里叶变换的代码实现如下所示（源代码的实现参考于这里）：
#define MY_PI 3.1415926535897932384626433
uint32_t bitReverse(uint32_t n, int dataArraySize)
{
uint32_t result = 0;
int processBitNum = log10(dataArraySize) / log10(2);
for (int i = 0; i < processBitNum; ++i)
{
if ((n & 1) == 1)
{
result = (result << 1) + 1;
}
else
{
result = result << 1;
}
n = n >> 1;
}
return result;
}
void myFFT(vector<float> & dataOfReal, vector<float> & dataOfImaginary)
{
int processStageNum = log10(dataOfReal.size()) / log10(2);
int processSampleNum = dataOfReal.size();
int exchangeIndex = 0;
float tempRealStart = 0.0f;
float tempImaginaryStart = 0.0f;
// Bit reverse order sorting
for (uint32_t i = 0; i < processSampleNum; i++)
{
exchangeIndex = bitReverse(i, processSampleNum);
if (i < exchangeIndex)
{
tempRealStart = dataOfReal[i];
dataOfReal[i] = dataOfReal[exchangeIndex];
dataOfReal[exchangeIndex] = tempRealStart;
tempImaginaryStart = dataOfImaginary[i];
dataOfImaginary[i] = dataOfImaginary[exchangeIndex];
dataOfImaginary[exchangeIndex] = tempImaginaryStart;
}
}
// Bit reverse order sorting
float realOfStageParameter;
float imaginaryOfStageParameter;
float tempReal;
float tempImaginary;
int blockStartNumOfStage;
int blockEndNumOfStage;
float realOfSinusoid;
float imaginaryOfSinusoid;
int jm1;
int ip;
for (int i = 1; i <= processStageNum; i++) //Loop for each processStageNum index layer
{
realOfStageParameter = 1.0f;
imaginaryOfStageParameter = 0.0f;
blockEndNumOfStage = pow(2, i);
blockStartNumOfStage = blockEndNumOfStage / 2;
realOfSinusoid = cos(MY_PI / blockStartNumOfStage);
imaginaryOfSinusoid = -sin(MY_PI / blockStartNumOfStage);
for (int j = 1; j <= blockStartNumOfStage; j++) //Loop for each sub DFT
{
jm1 = j - 1;
for (int k = jm1; k < processSampleNum; k = k + blockEndNumOfStage) //Loop for each butterfly
{
ip = k + blockStartNumOfStage;
tempReal = dataOfReal[ip] * realOfStageParameter - dataOfImaginary[ip] * imaginaryOfStageParameter;
tempImaginary = dataOfReal[ip] * imaginaryOfStageParameter + dataOfImaginary[ip] * realOfStageParameter;;
dataOfReal[ip] = dataOfReal[k] - tempReal;
dataOfImaginary[ip] = dataOfImaginary[k] - tempImaginary;
dataOfReal[k] = dataOfReal[k] + tempReal;
dataOfImaginary[k] = dataOfImaginary[k] + tempImaginary;
}
//update W parameter
tempReal = realOfStageParameter;
realOfStageParameter = tempReal * realOfSinusoid - imaginaryOfStageParameter * imaginaryOfSinusoid;
imaginaryOfStageParameter = tempReal * imaginaryOfSinusoid + imaginaryOfStageParameter * realOfSinusoid;
}
}
return;
}
由定义可知离散傅里叶变换和逆变换的定义如下所示：
- F [ n ] = ∑ k = 0 N f [ k ] e − j 2 π N k n F[n]=\sum\limits_{k=0}^{N}f[k]e^{{-j}{\frac{2\pi} {N}}kn} F[n]=k=0∑Nf[k]e−jN2πkn， n = 0 , 1 , 2 , . . , N − 1 n=0,1,2,..,N-1 n=0,1,2,..,N−1 ( D F T ) (DFT) (DFT)
- f [ k ] = 1 N ∑ n = 0 N F [ n ] e j 2 π N k n f[k]=\frac{1}{N}\sum\limits_{n=0}^{N}F[n]e^{{j}{\frac{2\pi} {N}}kn} f[k]=N1n=0∑NF[n]ejN2πkn， k = 0 , 1 , 2 , . . , N − 1 k=0,1,2,..,N-1 k=0,1,2,..,N−1 ( I D F T ) (IDFT) (IDFT)
如果现在已经求得了时间域信号的频率域信号 F [ n ] = a n + b n ∗ i F[n]=a_n+b_n*i F[n]=an+bn∗i，则时间域的信号为：
- f [ k ] = 1 N ∑ n = 0 N ( a n + b n ∗ j ) ∗ ( c o s ( X ) + j ∗ s i n ( X ) ) f[k]=\frac{1}{N}\sum\limits_{n=0}^{N}(a_n+b_n*j)*(cos(X)+j*sin(X)) f[k]=N1n=0∑N(an+bn∗j)∗(cos(X)+j∗sin(X))= 1 N ( a n ∗ c o s ( X ) − b n ∗ s i n ( X ) ) + 1 N ( a n ∗ s i n ( X ) + b n ∗ c o s ( X ) ) \frac{1}{N}(a_n*cos(X)-b_n*sin(X))+\frac{1}{N}(a_n*sin(X)+b_n*cos(X)) N1(an∗cos(X)−bn∗sin(X))+N1(an∗sin(X)+bn∗cos(X))
如果将求得的时间域信号的频率域信号 F [ n ] = a n + b n ∗ i F[n]=a_n+b_n*i F[n]=an+bn∗i的虚部取反 F [ n ] = a n − b n ∗ i F[n]=a_n-b_n*i F[n]=an−bn∗i然后以取反后的频率域信号 F [ n ] F[n] F[n]代替DFT公式中的 f [ k ] f[k] f[k]则可以得到：
- ∑ n = 0 N ( a n − b n ∗ j ) ∗ ( c o s ( X ) − j ∗ s i n ( X ) ) \sum\limits_{n=0}^{N}(a_n-b_n*j)*(cos(X)-j*sin(X)) n=0∑N(an−bn∗j)∗(cos(X)−j∗sin(X))= ( a n ∗ c o s ( X ) + b n ∗ s i n ( X ) ) − ( a n ∗ s i n ( X ) + b n ∗ c o s ( X ) ) (a_n*cos(X)+b_n*sin(X))-(a_n*sin(X)+b_n*cos(X)) (an∗cos(X)+bn∗sin(X))−(an∗sin(X)+bn∗cos(X))
从以上可以看出此时它们除了少一个乘积因子 1 N \frac{1}{N} N1之外就剩下虚部的符号相反而已，因此在求频率域信号 F [ n ] = a n + b n ∗ i F[n]=a_n+b_n*i F[n]=an+bn∗i的逆变换时可以先将其虚部取反 F [ n ] = a n − b n ∗ i F[n]=a_n-b_n*i F[n]=an−bn∗i然后带入计算DFT正变换的FFT代码最后将得到的结果的虚部取反并乘以乘积因子 1 N \frac{1}{N} N1就可以得到时间域信号。代码实现如下所示：
void myInverseFFT(vector<float>& dataOfReal, vector<float>& dataOfImaginary)
{
int processSampleNum = dataOfReal.size();
for (int i = 0; i < dataOfImaginary.size(); i++)
{
dataOfImaginary[i] = -dataOfImaginary[i];
}
myFFT(dataOfReal, dataOfImaginary);
for (int i = 0; i < dataOfImaginary.size(); i++)
{
dataOfReal[i] = dataOfReal[i]/ processSampleNum;
dataOfImaginary[i] = -dataOfImaginary[i]/ processSampleNum;
}
return;
}
以下是测试程序，测试程序的测试结果以及MATLAB对应的测试结果如图7和图8所示。
//测试程序
int main()
{
vector<float> dataOfReal;
vector<float> dataOfImaginary;
for (int i = 0; i < 8; i++)
{
dataOfReal.push_back(i + 1.0);
dataOfImaginary.push_back(0.0);
}
myFFT(dataOfReal, dataOfImaginary);
cout << "---------------------" << endl;
for (int i = 0; i < 8; i++)
{
cout << "real part=" << dataOfReal[i];
cout << " imaginary part=" << dataOfImaginary[i] << endl;
}
myInverseFFT(dataOfReal, dataOfImaginary);
cout << "---------------------" << endl;
for (int i = 0; i < 8; i++)
{
cout << "real part=" << dataOfReal[i];
cout << " imaginary part=" << dataOfImaginary[i] << endl;
}
}
如果直接按照离散傅里叶变换的公式直接计算，算法的时间复杂度为 O ( N 2 ) O(N^2) O(N2)，快速傅里叶变换的时间复杂度可以达到 O ( N ∗ log 2 N ) O(N*\log_2{N}) O(N∗log2N)。快速傅里叶变换算法有许多不同的变种。以上算法算是最常见的一种了，号称是Cooley–Tukey FFT algorithm。
接下来我们借助juce library来实现一个小程序。它读取输入通道的信号然后对该信号进行幅度谱的分析。该小程序是在这个库的一个demo：juce tutorial的基础上更改的，原本他使用的是其自带的快速傅里叶变换接口，现在我们使用我们自己的快速傅里叶变换接口来实现，顺便测试一下我们的FFT接口函数。我就不把全部代码放上来了，工程中下面两个函数中的IF语句实现了库的FFT函数和我们自己实现的FFT函数的切换。完整的工程代码下载地址为：https://download.csdn.net/download/caoleiwe/12673357
void drawNextFrameOfSpectrum()
{
if (0)//取0使用我们自身实现的FFT函数，取1使用库的FFT函数
{
// first apply a windowing function to our data
window.multiplyWithWindowingTable(fftData, fftSize);
// then render our FFT data..
forwardFFT.performFrequencyOnlyForwardTransform(fftData);
}
else
{
myFFT(dataOfReal, dataOfImaginary);
for (int i = 0; i < fftSize/2; i++)
{
fftData[i] = sqrt(dataOfReal[i]* dataOfReal[i] + dataOfImaginary[i] * dataOfImaginary[i]);
}
}
/
void pushNextSampleIntoFifo(float sample) noexcept
{
// if the fifo contains enough data, set a flag to say
// that the next frame should now be rendered..
if (fifoIndex == fftSize)
{
if (!nextFFTBlockReady)
{
zeromem(fftData, sizeof(fftData));
if (0) //取0使用我们自身实现的FFT函数，取1使用库的FFT函数
{
memcpy(fftData, fifo, sizeof(fifo));
}
else
{
for (int i = 0; i < fftSize; i++)
{
dataOfReal[i] = fifo[i];
dataOfImaginary[i] = 0.0f;
}
}
接下来我们来看一下效果。这里我们通过WG软件发送一个5000HZ的sine信号，然后查看其幅度谱的图形。FFT库函数的效果如图9所示。我们自己实现的FFT函数的效果如图10所示。图中右边为WG软件，左下角为WS软件接收的效果以及其对接收信号的幅度谱的分析。从图中可以看出我们自己实现的FFT函数的结果在除接收信号的其它频谱值的幅度比库函数分析的结果要高。这可能是因为我们只是简单的实现了FFT算法没有经过任何优化以及加窗函数（知识累计有限解释不了）等。
最后还要说一下这里频率的横坐标不是线性的，如果按线性的去画，对于48Khz的采样频率，6k到24k的范围就占了整个坐标范围的四分之三，这将极大的压缩低频信号的坐标范围。因此程序对频率坐标的范围进行了一个映射从而扩大了低频信号的坐标范围。工程中对应的语句如下：
auto skewedProportionX = 1.0f - std::exp(std::log(1.0f - i / (float)scopeSize) * 0.2f);
我们可以用MATLAB画两个图来理解。假设现在有索引 0 , 1 , 2 , 3 , 4 , 5 , . . . . , 500 0,1,2,3,4,5,....,500 0,1,2,3,4,5,....,500，第一个图以该索引为横坐标与纵坐标，第二个图的横坐标还是该索引但是纵坐标是该索引映射之后的值。结果如图11所示。我们可以看到映射之后0到300的索引范围（纵坐标）几乎映射到了百分之九十的横坐标的范围。
如果你觉得本文对你有帮助欢迎打赏。
1250
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
