1. 傅里叶变换说明
1.1. 中英文对照表
英文简称 | 英文全称 | 说明 |
---|---|---|
FT | Fourier Transformation | 傅里叶变换,针对 连续非周期 信号的傅里叶分析 |
FS | Fourier Serie | 傅里叶级数,任意 连续周期 信号在各正弦信号分量上的权重 |
DTFT | Discrete-time Fourier Transform | 离散时间傅里叶变换,离散非周期 信号的 FT |
DFT | Discrete Fourier Transform | 离散傅里叶变换,离散周期 信号的 FT |
FFT | Fast Fourier Transformation | 快速傅里叶变换 |
1.2 傅里叶相关变换公式
- FT (连续非周期信号):
F ( ω ) = ∫ − ∞ ∞ f ( t ) e − j ω t d t F(\omega)= \int_{-\infty}^{\infty}f(t)e^{-j{\omega}t}dt F(ω)=∫−∞∞f(t)e−jωtdt
- DTFT (连续非周期信号进行采样,然后做 FT ):
X ( ω ) = ∑ n = − ∞ + ∞ x ( n ) e − j ω n X(\omega) = \sum_{n=-\infty}^{+\infty}x(n)e^{-j\omega n} X(ω)=n=−∞∑+∞x(n)e−jωn
- DFT (在频域区间 [0, 2 π \pi π] 内以 2 π N \frac{2\pi}N N2π 为间隔对 DTFT 进行采样):
X ( k ) = ∑ n = 0 N − 1 x ( n ) e − j 2 π N k n X(k) = \sum_{n=0}^{N-1}x(n)e^{-j\frac{2\pi}Nkn} X(k)=n=0∑N−1x(n)e−jN2πkn
-
IDFT (DFT 反变换,注意系数 1 N \frac1N N1)
x ( n ) = 1 N ∑ k = 0 N − 1 X ( k ) e j 2 π N n k x(n) = \frac1N\sum_{k=0}^{N-1}X(k)e^{j\frac{2\pi}{N}nk} x(n)=N1k=0∑N−1X(k)ejN2πnk -
FFT
-
相关性质
令 W N k n = e − j 2 π N k n W_N^{kn} = e^{-j\frac{2\pi}Nkn} WNkn=e−jN2πkn , 则:
1). W N 0 = W N N = 1 W_N^0 = W_N^N = 1 WN0=WNN=1
2). W N k = W N k + N W_N^k = W_N^{k+N} WNk=WNk+N
3). W d N d k = W N k W_{dN}^{dk} = W_{N}^{k} WdNdk=WNk
4). W 2 1 = − 1 W_2^1 = -1 W21=−1
5). W N k + N 2 = − W N k W_N^{k+\frac{N}2} = -W_N^k WNk+2N=−WNk
6). W N k + n = W N k ∗ W N n W_N^{k+n} = W_N^k * W_N^n WNk+n=WNk∗WNn
-
公式推导
设有如下多项式,其中 n n n 为 2 的整数次幂:
f ( x ) = a 0 + a 1 x 1 + a 2 x 2 + a 3 x 3 + . . . + a n − 1 x n − 1 f(x) = a_0 + a_1x^1 + a_2x^2 + a_3x^3 + ... + a_{n-1}x^{n-1} f(x)=a0+a1x1+a2x2+a3x3+...+an−1xn−1
分离偶数次幂项与奇数次幂项,并令:
f 1 ( x ) = a 0 + a 2 x 2 + a 4 x 4 + . . . + a n − 2 x n − 2 f_1(x) = a_0 + a_2x^2 + a_4x^4 + ... + a_{n-2}x^{n-2} f1(x)=a0+a2x2+a4x4+...+an−2xn−2f 2 ( x ) = a 1 + a 3 x 2 + a 5 x 4 + . . . + a n − 1 x n − 2 f_2(x) = a_1 + a_3x^2 + a_5x^4 + ... + a_{n-1}x^{n-2} f2(x)=a1+a3x2+a5x4+...+an−1xn−2
则有:
f ( x ) = f ′ ( x 2 ) = f 1 ( x 2 ) + x f 2 ( x 2 ) f(x) = f'(x^2)= f_1(x^2) + xf_2(x^2) f(x)=f′(x2)=f1(x2)+xf2(x2)
令 x = w n k x = w_n^k x=wnk ,且 k < n 2 k < \frac{n}2 k<2n 时,由性质 3 可得
f ( w n k ) = f ′ ( ( w n k ) 2 ) = f 1 ( ( w n k ) 2 ) + w n k f 2 ( ( w n k ) 2 ) = f 1 ( ( w n / 2 k ) ) + w n k f 2 ( ( w n / 2 k ) ) f(w_n^k) = f'((w_n^k)^2) = f_1((w_n^k)^2) + w_n^kf_2((w_n^k)^2) = f_1((w_{n/2}^k)) + w_n^kf_2((w_{n/2}^k)) f(wnk)=f′((wnk)2)=f1((wnk)2)+wnkf2((wnk)2)=f1((wn/2k))+wnkf2((wn/2k))
当 k ′ k' k′ 取后半部分,即 k ′ > = n 2 k' >= \frac{n}2 k′>=2n 时,由性质 2 、 5 得
f ( w n k ′ ) = f 1 ( ( w n k + n / 2 ) 2 ) + w n k + n / 2 f 2 ( ( w n k + n / 2 ) 2 ) = f 1 ( ( w n / 2 k ) ) − w n k f 2 ( ( w n / 2 k ) ) f(w_n^{k'}) = f_1((w_n^{k+{n/2}})^2) + w_n^{k+{n/2}}f_2((w_n^{k+{n/2}})^2) = f_1((w_{n/2}^k)) - w_n^kf_2((w_{n/2}^k)) f(wnk′)=f1((wnk+n/2)2)+wnk+n/2f2((wnk+n/2)2)=f1((wn/2k))−wnkf2((wn/2k))
结合 DFT 公式,可知1)n 点的 DFT 可以转换成 n/2 的 DFT
2)计算前一半 DFT 值时,即可同时得到后一半 DFT 值
-
矩阵示例
对信号 x ( t ) x(t) x(t) 在一个周期内进行 N N N 次采样,则 N N N 点的 DFT 满足:
{ w n 00 w n 01 ⋯ w n 0 ( n − 1 ) w n 10 w n 11 ⋯ w n 1 ( n − 1 ) ⋮ ⋮ ⋱ ⋮ w n ( n − 1 ) 0 w n ( n − 1 ) 1 ⋯ w n ( n − 1 ) ( n − 1 ) } ∗ { x 0 x 1 ⋮ x n − 1 } = { X 0 X 1 ⋮ X n − 1 } \begin{Bmatrix} w_n^{00}&w_n^{01}&\cdots&w_n^{0(n-1)}\\ w_n^{10}&w_n^{11}&\cdots&w_n^{1(n-1)}\\ \vdots&\vdots&\ddots&\vdots\\ w_n^{(n-1)0}&w_n^{(n-1)1}&\cdots&w_n^{(n-1)(n-1)} \end{Bmatrix} * \begin{Bmatrix} x_0\\ x_1\\ \vdots\\ x_{n-1} \end{Bmatrix} = \begin{Bmatrix} X_0\\ X_1\\ \vdots\\ X_{n-1} \end{Bmatrix} ⎩ ⎨ ⎧wn00wn10⋮wn(n−1)0wn01wn11⋮wn(n−1)1⋯⋯⋱⋯wn0(n−1)wn1(n−1)⋮wn(n−1)(n−1)⎭ ⎬ ⎫∗⎩ ⎨ ⎧x0x1⋮xn−1⎭ ⎬ ⎫=⎩ ⎨ ⎧X0X1⋮Xn−1⎭ ⎬ ⎫
n n n 表示采样点下标,即 w w w 矩阵列数; k k k 表示周期内的频率分量下标,即 w w w 矩阵行数
将 w w w 矩阵偶、奇列重排,DFT 等式可表示为
W ′ ∗ { x 0 x 2 ⋮ x n − 2 x 1 x 3 ⋮ x n − 1 } = { X 0 X 2 ⋮ X n − 2 X 1 X 3 ⋮ X n − 1 } W'* \begin{Bmatrix} x_0\\ x_2\\ \vdots\\ x_{n-2}\\ x_1\\ x_3\\ \vdots\\ x_{n-1} \end{Bmatrix} = \begin{Bmatrix} X_0\\ X_2\\ \vdots\\ X_{n-2}\\ X_1\\ X_3\\ \vdots\\ X_{n-1} \end{Bmatrix} W′∗⎩ ⎨ ⎧x0x2⋮xn−2x1x3⋮xn−1⎭ ⎬ ⎫=⎩ ⎨ ⎧X0X2⋮Xn−2X1X3⋮Xn−1⎭ ⎬ ⎫
此时,
W ′ = { w n 00 w n 02 ⋯ w n 0 ( n − 2 ) w n 01 w n 03 ⋯ w n 0 ( n − 1 ) w n 10 w n 12 ⋯ w n 1 ( n − 2 ) w n 11 w n 13 ⋯ w n 1 ( n − 1 ) ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ ⋱ ⋮ w n ( n 2 − 1 ) 0 w n ( n 2 − 1 ) 2 ⋯ w n ( n 2 − 1 ) ( n − 2 ) w n ( n 2 − 1 ) 1 w n ( n 2 − 1 ) 3 ⋯ w n ( n 2 − 1 ) ( n − 1 ) w n n 2 0 w n n 2 2 ⋯ w n n 2 ( n − 2 ) w n n 2 1 w n n 2 3 ⋯ w n n 2 ( n − 1 ) w n ( n 2 + 1 ) 0 w n ( n 2 + 1 ) 2 ⋯ w n ( n 2 + 1 ) ( n − 2 ) w n ( n 2 + 1 ) 1 w n ( n 2 + 1 ) 3 ⋯ w n ( n 2 + 1 ) ( n − 1 ) ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ ⋱ ⋮ w n ( n − 1 ) 0 w n ( n − 1 ) 2 ⋯ w n ( n − 1 ) ( n − 2 ) w n ( n − 1 ) 1 w n ( n − 1 ) 3 ⋯ w n ( n − 1 ) ( n − 1 ) } W' = \begin{Bmatrix} w_n^{00}&w_n^{02}&\cdots&w_n^{0(n-2)}&w_n^{01}&w_n^{03}&\cdots&w_n^{0(n-1)}\\ w_n^{10}&w_n^{12}&\cdots&w_n^{1(n-2)}&w_n^{11}&w_n^{13}&\cdots&w_n^{1(n-1)}\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\ w_n^{(\frac{n}2-1)0}&w_n^{(\frac{n}2-1)2}&\cdots&w_n^{(\frac{n}2-1)(n-2)}&w_n^{(\frac{n}2-1)1}&w_n^{(\frac{n}2-1)3}&\cdots&w_n^{(\frac{n}2-1)(n-1)}\\ w_n^{\frac{n}20}&w_n^{\frac{n}22}&\cdots&w_n^{\frac{n}2(n-2)}&w_n^{\frac{n}21}&w_n^{\frac{n}23}&\cdots&w_n^{\frac{n}2(n-1)}\\ w_n^{(\frac{n}2+1)0}&w_n^{(\frac{n}2+1)2}&\cdots&w_n^{(\frac{n}2+1)(n-2)}&w_n^{(\frac{n}2+1)1}&w_n^{(\frac{n}2+1)3}&\cdots&w_n^{(\frac{n}2+1)(n-1)}\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\ w_n^{(n-1)0}&w_n^{(n-1)2}&\cdots&w_n^{(n-1)(n-2)}&w_n^{(n-1)1}&w_n^{(n-1)3}&\cdots&w_n^{(n-1)(n-1)} \end{Bmatrix} W′=⎩ ⎨ ⎧wn00wn10⋮wn(2n−1)0wn2n0wn(2n+1)0⋮wn(n−1)0wn02wn12⋮wn(2n−1)2wn2n2wn(2n+1)2⋮wn(n−1)2⋯⋯⋱⋯⋯⋯⋱⋯wn0(n−2)wn1(n−2)⋮wn(2n−1)(n−2)wn2n(n−2)wn(2n+1)(n−2)⋮wn(n−1)(n−2)wn01wn11⋮wn(2n−1)1wn2n1wn(2n+1)1⋮wn(n−1)1wn03wn13⋮wn(2n−1)3wn2n3wn(2n+1)3⋮wn(n−1)3⋯⋯⋱⋯⋯⋯⋱⋯wn0(n−1)wn1(n−1)⋮wn(2n−1)(n−1)wn2n(n−1)wn(2n+1)(n−1)⋮wn(n−1)(n−1)⎭ ⎬ ⎫将 W ′ W' W′ 矩阵划分为四部分
W ′ = { A B C D } W' = \begin{Bmatrix} A&B\\ C&D \end{Bmatrix} W′={ACBD}其中, A A A 矩阵中各项为
A = { w n 00 w n 02 ⋯ w n 0 ( n − 2 ) w n 10 w n 12 ⋯ w n 1 ( n − 2 ) ⋮ ⋮ ⋱ ⋮ w n ( n 2 − 1 ) 0 w n ( n 2 − 1 ) 2 ⋯ w n ( n 2 − 1 ) ( n − 2 ) } A = \begin{Bmatrix} w_n^{00}&w_n^{02}&\cdots&w_n^{0(n-2)}\\ w_n^{10}&w_n^{12}&\cdots&w_n^{1(n-2)}\\ \vdots&\vdots&\ddots&\vdots\\ w_n^{(\frac{n}2-1)0}&w_n^{(\frac{n}2-1)2}&\cdots&w_n^{(\frac{n}2-1)(n-2)} \end{Bmatrix} A=⎩ ⎨ ⎧wn00wn10⋮wn(2n−1)0wn02wn12⋮wn(2n−1)2⋯⋯⋱⋯wn0(n−2)wn1(n−2)⋮wn(2n−1)(n−2)⎭ ⎬ ⎫则 1)依据性质 2 可得
C = A C = A C=A2)依据性质 6 可得(注意此处 k < n / 2 k < n/2 k<n/2 )
B = { w n 0 ⋯ ⋯ ⋯ ⋯ w n 1 ⋯ ⋯ ⋯ ⋯ w n k ⋯ ⋯ ⋯ ⋯ w n ( n 2 − 1 ) } ∗ A = d i a g ( w n k ) ∗ A B = \begin{Bmatrix} w_n^{0}&\cdots&\cdots&\cdots\\ \cdots&w_n^{1}&\cdots&\cdots\\ \cdots&\cdots&w_n^{k}&\cdots\\ \cdots&\cdots&\cdots&w_n^{(\frac{n}2-1)}\\ \end{Bmatrix} *A = diag(w_n^k)*A B=⎩ ⎨ ⎧wn0⋯⋯⋯⋯wn1⋯⋯⋯⋯wnk⋯⋯⋯⋯wn(2n−1)⎭ ⎬ ⎫∗A=diag(wnk)∗A
3)依据性质 2 5 可得
D = − B = − d i a g ( w n k ) ∗ A D = -B = -diag(w_n^k)*A D=−B=−diag(wnk)∗A4)依据性质 3 可将 A A A 矩阵简化为
A = { w n 2 00 w n 2 01 ⋯ w n 2 0 ( n 2 − 1 ) w n 2 10 w n 2 11 ⋯ w n 2 1 ( n 2 − 1 ) ⋮ ⋮ ⋱ ⋮ w n 2 ( n 2 − 1 ) 0 w n 2 ( n 2 − 1 ) 1 ⋯ w n 2 ( n 2 − 1 ) ( n 2 − 1 ) } A = \begin{Bmatrix} w_{\frac{n}2}^{00}&w_{\frac{n}2}^{01}&\cdots&w_{\frac{n}2}^{0(\frac{n}2-1)}\\ w_{\frac{n}2}^{10}&w_{\frac{n}2}^{11}&\cdots&w_{\frac{n}2}^{1(\frac{n}2-1)}\\ \vdots&\vdots&\ddots&\vdots\\ w_{\frac{n}2}^{(\frac{n}2-1)0}&w_{\frac{n}2}^{(\frac{n}2-1)1}&\cdots&w_{\frac{n}2}^{(\frac{n}2-1)(\frac{n}2-1)} \end{Bmatrix} A=⎩ ⎨ ⎧w2n00w2n10⋮w2n(2n−1)0w2n01w2n11⋮w2n(2n−1)1⋯⋯⋱⋯w2n0(2n−1)w2n1(2n−1)⋮w2n(2n−1)(2n−1)⎭ ⎬ ⎫
此时 A A A 即 n / 2 n/2 n/2 点 DFT 的 W n 2 W_{\frac{n}2} W2n ,对于 n n n 点的 DFT 可以拆分为 n / 2 n/2 n/2 点的 DFT ,且系数因子(即旋转因子)满足:
W n = { W n 2 d i a g ( w n k ) ∗ W n 2 W n 2 − d i a g ( w n k ) ∗ W n 2 } W_n = \begin{Bmatrix} W_{\frac{n}2}&diag(w_n^k)*W_{\frac{n}2}\\ W_{\frac{n}2}&-diag(w_n^k)*W_{\frac{n}2} \end{Bmatrix} Wn={W2nW2ndiag(wnk)∗W2n−diag(wnk)∗W2n}特别地,对于 2 点的 DFT 有
{ X 0 X 1 } = W 2 ∗ { x 0 x 1 } = { 1 1 1 − 1 } ∗ { x 0 x 1 } \begin{Bmatrix} X_0\\ X_1 \end{Bmatrix} = W_2 * \begin{Bmatrix} x_0\\ x_1 \end{Bmatrix} = \begin{Bmatrix} 1&1\\ 1&-1 \end{Bmatrix} * \begin{Bmatrix} x_0\\ x_1 \end{Bmatrix} {X0X1}=W2∗{x0x1}={111−1}∗{x0x1}8 点 DFT 蝶形运算图示
-
代码实现
function y=myditfft (x) %本程序对输入序列 x实现DIT-FFT基2算法,点数取大于等于x长度的2的幂次 %y=myditfft (x) x=randn(1,2^4); m=nextpow2 (length(x)); %进行补零操作,满足2的整数幂----nextpow2(7)=3, N=2 ^m; %求x的长度对应的2的最低幂次m,对应多少点的FFT if length(x) <N x=[x, zeros(1, N-length(x))] ;%若x的长度不是2的幂,补零到2的整数幂. end nxd=bin2dec (fliplr (dec2bin((1:N)-1,m)))+1; %求1:2^m数列的倒序 %bin2dec将二进制转换为10进制 flip1r倒叙排列 dec2bin把一个十进制数转换成一个字符串形式的二进制数 y=x(nxd) ; %将x倒序排列作为y的初始值,这样我们输出的数值和蝶形运算的数值相对应。 for mm=1 :m %将DFT作m次基2分解,从左到右,对每次分解作DFT运算——把DFT分成log2N段来实现 Nmr= 2 ^mm; %这里是旋转因子中的N u=1; %初始化WN0=1 WN=exp (-1i*2*pi/Nmr) ; %本次分解的基本DFT因子WN=exp (-i*2*pi/ Nmr) for j=1 :Nmr/2 %本次跨越间隔内的各次蝶形运算 for k=j :Nmr :N %本次蝶形运算的跨越间隔为Nmr=2 ^mm kp=k+Nmr/2; %确定蝶形运算的对应单元下标. t=y (kp) *u; %蝶形运算的乘积项 y (kp)=y(k) -t; %蝶形运算的减法项 y(k)=y(k) +t; %蝶形运算的加法项 end u=u*WN; %修改旋转因子,多乘一个基本DFT因子WN end end
【参考1】: FFT 代码实现
【参考2】: 蝶形运算输入数据重排原理
【参考3】: 图示蝶形运算
-
1.3 其他图示
-
矩形波与正弦波
-
公式
y ( x ) = ∑ n = 0 N 1 2 n + 1 s i n ( ( 2 n + 1 ) x ) y(x) = \sum_{n=0}^N \frac1{2n+1}sin((2n+1)x) y(x)=n=0∑N2n+11sin((2n+1)x) -
图示
-