离散分数傅立叶转换

离散分数傅立叶转换（Discrete Fractional Fourier Transform）是用来解决数字序列分数傅立叶转换的计算问题，方法是利用它们的特征函数展开的表达来实现离散算法，而离散分数傅立叶转换的特征函数是埃尔米特多项式与高斯函数的乘积，这样的特征函数同时也是傅立叶转换的特征函数。利用离散傅立叶转换的结果，可以建立周期分数傅立叶转换的离散算法。

## 离散傅立叶转换及其周期性质

在建立离散分数傅立叶转换的算法之前，必须对离散傅立叶转换（Discrete Fourier Transform）有一定的理论认识。

对于一般的傅立叶转换具有4周期性质，也就是对任何讯号连续做4次傅立叶转换就会回到它自己。然而离散傅立叶转换也有这样的性质，所以以下先讨论离散傅立叶转换的4周期性质。

### 离散傅立叶转换的矩阵

对于一串数字序列讯号$x_0,x_1,...,x_{N-1}$，定义它的离散傅立叶转换是：

$$y_m=\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1} x_ne^{-\frac{2\pi mni}{N}} ,m=0,1,...,N-1$$

引入矩阵符号

$$X=\begin{bmatrix} x_0 \\ x_1 \\ \vdots \\ x_{N-1} \end{bmatrix} \quad Y=\begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_{N-1} \end{bmatrix}$$

$$E=\frac{1}{\sqrt{N}}\begin{bmatrix} W^{0 \times 0} & W^{0 \times 1} & \vdots & W^{0 \times N-1} \\ W^{1 \times 0} & W^{1 \times 1} & \vdots & W^{1 \times N-1} \\ \vdots & \vdots & \ddots & \vdots \\ W^{N-1 \times 0} & W^{N-1 \times 1} & \vdots & W^{N-1 \times N-1} \\ \end{bmatrix}$$

其中 $W=e^{-\frac{\pi 2i}{\sqrt{N}}}$，现在可以把离散傅立叶转换写成： $Y=EX$

由上面的分析可知一串长度为N的数字序列讯号$x_0,x_1,...,x_{N-1}$，其离散傅立叶转换在N维向量空间相当于一个线性转换，也就是矩阵E，矩阵E称作离散傅立叶转换的矩阵。又因为

$$\frac{1}{N}\sum_{l=0}^{N-1}W^{l \times m} \times \overline{W}^{l \times n}=\frac{1}{N}\sum_{l=0}^{N-1}e^{-\frac{2\pi li}{N} (m-n)}= \begin{cases} 1 \quad m=n \\ 0 \quad m \ne n \end{cases}$$

所以 $EE^*=I$

离散傅立叶转换的线性转换是正交转换。

### 离散傅立叶转换4周期性质的证明

利用离散傅立叶转换的矩阵，可以证明4周期性质。因为

$$\frac{1}{N}\sum_{l=0}^{N-1}W^{l \times m} \times \overline{W}^{l \times n}=\frac{1}{N}\sum_{l=0}^{N-1}e^{-\frac{2\pi li}{N} (m-n)}= \begin{cases} 1 \quad \mod(m+n)=0 \\ 0 \quad \mod(m+n) \ne 0 \end{cases}$$

而且$0\leqslant m,n \leqslant N-1$，所以

$$E^2=\begin{bmatrix} 1 & 0 & \vdots & 0 \\ 0 & 0 & \vdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 1 & \vdots & 0 \end{bmatrix}$$

由上可知讯号$x_0,x_1,...,x_{N-1}$做两次离散傅立叶转换后会得到$x_0,x_{N-1},...,x_1$，除了第一个数字不动其余顺序皆颠倒，这意味着发生了时间反射。现在用矩阵来表达连续两次的离散傅立叶转换：

$$Z=E^2X=\begin{bmatrix} 1 & 0 & \vdots & 0 \\ 0 & 0 & \vdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 1 & \vdots & 0 \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \\ \vdots \\ x_{N-1} \end{bmatrix}=\begin{bmatrix} x_0 \\ x_{N-1} \\ x_{N-2} \\ \vdots \\ x_1 \\ \end{bmatrix}$$

再来计算连续三次的转换：

$$E^3=E^2E=\frac{1}{\sqrt{N}}\begin{bmatrix} 1 & 0 & \vdots & 0 \\ 0 & 0 & \vdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 1 & \vdots & 0 \end{bmatrix} \begin{bmatrix} W^{0 \times 0} & W^{0 \times 1} & \vdots & W^{0 \times N-1} \\ W^{1 \times 0} & W^{1 \times 1} & \vdots & W^{1 \times N-1} \\ \vdots & \vdots & \ddots & \vdots \\ W^{N-1 \times 0} & W^{N-1 \times 1} & \vdots & W^{N-1 \times N-1} \\ \end{bmatrix}$$

$$=\frac{1}{\sqrt{N}}\begin{bmatrix} W^{0 \times 0} & W^{0 \times 1} & \vdots & W^{0 \times N-1} \\ W^{N-1 \times 0} & W^{N-1 \times 1} & \vdots & W^{N-1 \times N-1} \\ \vdots & \vdots & \ddots & \vdots \\ W^{1 \times 0} & W^{1 \times 1} & \vdots & W^{1 \times N-1} \\ \end{bmatrix}$$

最后计算连续四次的转换：

$$E^4=E^2 \times E^2 = \begin{bmatrix} 1 & 0 & \vdots & 0 \\ 0 & 0 & \vdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 1 & \vdots & 0 \end{bmatrix} \begin{bmatrix} 1 & 0 & \vdots & 0 \\ 0 & 0 & \vdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 1 & \vdots & 0 \end{bmatrix}=\begin{bmatrix} 1 & 0 & \vdots & 0 \\ 0 & 1 & \vdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \vdots & 1 \end{bmatrix}$$

由上式可知，离散傅立叶转换连续做4次一定会让讯号返回它自己，这就是离散傅立叶转换4周期性质的证明。

## 离散傅立叶转换的周期性定理

为了得到离散傅立叶转换的周期性定理，必须理解它的运算性质，因此现在用另外一个观点去说明4周期性质。

$$U=E^3X=E \times E^2X$$

$$=\begin{bmatrix} W^{0 \times 0} & W^{0 \times 1} & \vdots & W^{0 \times N-1} \\ W^{1 \times 0} & W^{1 \times 1} & \vdots & W^{1 \times N-1} \\ \vdots & \vdots & \ddots & \vdots \\ W^{N-1 \times 0} & W^{N-1 \times 1} & \vdots & W^{N-1 \times N-1} \\ \end{bmatrix} \begin{bmatrix} 1 & 0 & \vdots & 0 \\ 0 & 0 & \vdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 1 & \vdots & 0 \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ \vdots \\ x_{N-1} \end{bmatrix}$$

$$=\begin{bmatrix} W^{0 \times 0} & W^{0 \times 1} & \vdots & W^{0 \times N-1} \\ W^{1 \times 0} & W^{1 \times 1} & \vdots & W^{1 \times N-1} \\ \vdots & \vdots & \ddots & \vdots \\ W^{N-1 \times 0} & W^{N-1 \times 1} & \vdots & W^{N-1 \times N-1} \\ \end{bmatrix} \begin{bmatrix} x_0 \\ x_{N-1} \\ x_{N-2} \\ \vdots \\ x_1 \\ \end{bmatrix}$$

$$=\begin{bmatrix} 1 & 0 & \vdots & 0 \\ 0 & 0 & \vdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 1 & \vdots & 0 \end{bmatrix} \begin{bmatrix} W^{0 \times 0} & W^{0 \times 1} & \vdots & W^{0 \times N-1} \\ W^{1 \times 0} & W^{1 \times 1} & \vdots & W^{1 \times N-1} \\ \vdots & \vdots & \ddots & \vdots \\ W^{N-1 \times 0} & W^{N-1 \times 1} & \vdots & W^{N-1 \times N-1} \\ \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ \vdots \\ x_{N-1} \end{bmatrix}$$

$$=\begin{bmatrix} y_0 \\ y_{N-1} \\ y_{N-2} \\ \vdots \\ y_1 \end{bmatrix}$$

上面的分析可以理解为原始讯号作三次离散傅立叶转换其成分还是跟原本讯号相同，只是顺序上做了调整，也就是从第二个讯号至最后一个讯号作了颠倒的顺序。又或者可以这样理解，原本的讯号经过交换顺序而得到$x_0,x_{N-1},x_{N-2},...,x_1$，其离散傅立叶转换就是$y_0,y_{N-1},y_{N-2},...,y_1$，若再做一次离散傅立叶转换，结果就是将$x_0,x_{N-1},x_{N-2},...,x_1$的下标按照前面规则作了交换，连续四次的离散傅立叶转换其结果应该为$x_0,x_1,...,x_{N-1}$，由此可知离散傅立叶转换具有4周期的性质。

因此，$E^4=1$，即任何数字序列作四次离散傅立叶转换就会返回它自己，在矩阵乘法之下，集合 $\begin{Bmatrix} I,E,E^2,E^3 \end{Bmatrix}$ 会构成一个四阶循环群。

### 定理

离散傅立叶转换的整数次重复运算等价于一个4阶循环矩阵群 $\begin{Bmatrix} I,E,E^2,E^3 \end{Bmatrix}$。

## 离散分数傅立叶转换的计算方法

利用前述离散傅立叶转换的计算方法可以得到离散分数傅立叶转换的计算方法，因此接下来的分析就先从离散傅立叶转换出发。

### 推导

根据离散傅立叶转换的矩阵形式，定义符号：

$$X^{(l)}=\begin{bmatrix} x_0^{(l)} \\ x_1^{(l)} \\ \vdots \\ x_{N-1}^{(l)} \end{bmatrix} ,l=0,1,2,3$$

其递回关系为：

$$X^{(l+1)}=EX^{(l)}=E^{l+1}X^{(0)},l=0,1,2,3$$

当 $l=0$ 时，对应的是原始数字序列讯号；当 $l=1$ 时，对应的昰原始数字序列讯号的离散傅立叶转换。现在定义幂次a的离散分数傅立叶转换为：

$$X^{(a)}=\sum_{l=0}^{3}A_l(a)X^{(l)}=(\sum_{l=0}^{3}A_l(a)E^l)X^{(0)}$$

上式中的 $A_j(a),j=0,1,2,3$ 是a的连续函数。引入矩阵符号为：

$$E^a=\sum_{l=0}^{3}A_l(a)E^l$$

称作矩阵E的a次幂，因此幂次a的离散分数傅立叶转换可以写成：

$$X^{(a)}=E^aX^{(0)}$$

其中 $A_j(a),j=0,1,2,3$ 是只与a有关的系数，它使得 $a \in R^+,F_s^a$ 全体生成的矩阵族 $\begin{Bmatrix} F_s^a ,a \in R^+ \end{Bmatrix}$ 满足以下公理。

- 连续性公理： $\forall a \in R,E^a:C^N \to C^N$ 为连续的。
- 边界性公理： $\forall a \in R,E^a$ 是离散傅立叶转换连续作用a次的结果。
- 周期性公理： $\forall a \in R,n \in N,E^a=E^{a+4n}$。
- 可加性公理： $\forall a,b \in R,E^a E^b = E^b E^a = E^{a+b}$。

由以上可得连续函数 $A_j(a),A_i(b),i,j = 0,1,2,3$ 的方程组：

$$\begin{cases} A_0(a+b)=A_0(a)A_0(b)+A_1(a)A_3(b)+A_2(a)A_2(b)+A_3(a)A_1(b) \\ A_1(a+b)=A_0(a)A_1(b)+A_1(a)A_0(b)+A_2(a)A_3(b)+A_3(a)A_2(b) \\ A_2(a+b)=A_0(a)A_2(b)+A_1(a)A_1(b)+A_2(a)A_0(b)+A_3(a)A_3(b) \\ A_3(a+b)=A_0(a)A_3(b)+A_1(a)A_2(b)+A_2(a)A_1(b)+A_3(a)A_0(b) \end{cases}$$

为了解得这个方程组的解，先做变数变换：

$$\begin{bmatrix} B_0(a) \\ B_1(a) \\ B_2(a) \\ B_3(a) \end{bmatrix}=\begin{bmatrix} 1 & i & 1 & 1 \\ 1 & -i & -1 & i \\ 1 & 1 & 1 & -1 \\ 1 & i & -1 & -i \end{bmatrix} \begin{bmatrix} A_0(a) \\ A_1(a) \\ A_2(a) \\ A_3(a) \end{bmatrix}$$

则方程组可表达为： $B_l(a+b)=B_l(a)B_l(b),l=0,1,2,3$

同时又由边界性公理可得： $A_(m)=\delta(m-l)= \begin{cases} 1 \quad m=l \\ 0 \quad m \ne l \end{cases}$

因此 $B_l(a),l=0,1,2,3$ 的边界条件就是： $B_l(m)=e^{-\frac{2\pi lmi}{4}}$

此外周期性条件为： $B_l(a+4)=B_l(a)$

考虑了这些条件后方程组 $B_l(a+b)=B_l(a)B_l(b),l=0,1,2,3$ 有唯一解：

$$B_l(a)=e^{-\frac{2\pi lai}{4}},j=0,1,2,3$$

最后再根据

$$\begin{bmatrix} B_0(a) \\ B_1(a) \\ B_2(a) \\ B_3(a) \end{bmatrix}=\begin{bmatrix} 1 & i & 1 & 1 \\ 1 & -i & -1 & i \\ 1 & 1 & 1 & -1 \\ 1 & i & -1 & -i \end{bmatrix} \begin{bmatrix} A_0(a) \\ A_1(a) \\ A_2(a) \\ A_3(a) \end{bmatrix}$$

得到方程组的通解为：

$$A_l(a)=\frac{1}{4} \frac{1-e^{-2\pi (a-l)i}}{1-e^{-\frac{2\pi (a-l)i}{4}}}=\cos(\frac{(a-l)\pi}{4})\cos(\frac{2(a-l)\pi}{4})e^{-\frac{3i(a-l)\pi}{4}},l=0,1,2,3$$

所以满足前述4个公理的矩阵E之a次幂定义为：

$$E^a=\sum_{l=0}^{3}\cos(\frac{(a-l)\pi}{4})\cos(\frac{2(a-l)\pi}{4})e^{-\frac{3i(a-l)\pi}{4}}E^l$$

### 结果

上式得到了 $E^a$ 的离散形式并把它带入 $X^{(a)}=E^aX^{(0)}$ 最后得到了离散分数傅立叶转换的计算方法：

$$X^{(a)}=(\sum_{l=0}^{3}\cos(\frac{(a-l)\pi}{4})\cos(\frac{2(a-l)\pi}{4})e^{-\frac{3i(a-l)\pi}{4}}E^l)X^{(0)}$$

若去对照分数傅立叶转换的形式，上式的计算方法可以说是C.C.Shih的分数傅立叶转换的离散算法。