離散分數傅立葉轉換

離散分數傅立葉轉換（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的分數傅立葉轉換的離散演算法。