一點說明
1. 完成情況
- 傅里葉變換（DFT）與余弦變換（DCT）
- 快速傅里葉變換（FFT）與快速傅里葉逆變換（IFFT）
- 快速余弦變換（FCT）與快速余弦逆變換（IFCT）
2. 項目結構
- root/
- bin/（存放可執行文件，雙擊直接運行）
- src/（存放項目源代碼，編譯需要先配置OpenCV）
- Photos/（存放圖片樣例，包括原圖與變換后的圖片）
- Data/（存放圖像變換生成的數據，如幅度值與相位值）
- README.md（技術分析報告）
思路分析
1. 傅里葉變換
利用二維的傅里葉變換，我們可以將圖像信號從空間域（spatial domain）變換到其對應的頻域（frequency domain）中進行分析，這與我們日常觀察世界的視角是截然不同的。
所謂傅里葉變換其實是正交變換的一種，其原理是“周期與非周期信號都可用正弦函數的加權積分表示”。在圖像處理中，我們一般用到的是二維離散傅里葉變換，具體公式如下。
在具體的編碼實現時，我們常常利用歐拉變換將公式中的實部與虛部分離（如下圖）。
可以發現， 通過上述變換，我們將二維離散傅里葉變換的公式由\(F(u,v)=|F(u,v)|e^{jφ(u,v)}\)的形式轉化成了\(F(u,v)=R(u,v)+jI(u,v)\)的形式，所以有下式。
\[φ(u,v)=arctan[\frac{I(u,v)}{R(u,v)}]\\ |F(u,v)|=[R^2(u,v)+I^2(u,v)]^{\frac{1}{2}} \]
上式中的 \(φ(u,v)\)為圖像信號在\((u,v)\)點的幅度值，\(|F(u,v)|\)為信號在\((u,v)\)點的相位值。通過上面的拆分，我們可以輕松地編寫程序進行計算，得到所謂的幅度圖像和相位圖像。其中幅度圖像包含了我們所需要的圖像幾何結構的所有信息，在圖像分析和處理中應用最為廣泛。
2. 余弦變換
余弦變換也是正交變換的一種。由前面提到的傅里葉變換公式可以知道，偶函數的傅立葉變換的虛部為零，因而不需要計算，只計算余弦項變換，這就是余弦變換 。顯然，余弦變換的變換核為實數的余弦函數（公式如下），其計算速度相比傅里葉變換要快得多，所以余弦變換被廣泛應用於圖像有損壓縮和語音信號處理等眾多領域。
3. 逆變換
可以證明，傅里葉變換與余弦變換都是可逆的，其逆變換公式可以通過矩陣的逆運算性質求出，具體公式如下。
4. 快速變換
無論是快速傅里葉變換還是快速余弦變換，其實都是利用了變換的可分離性，借助動態規划的思想，將二維離散變換的復雜度從 \(O(n^2)\) 優化到 \(O(nlogn)\) 。基本思路如下圖所示（圖源博客）。
關鍵代碼
此處僅貼出了快速傅里葉正逆變換(FFT&IFFT)與快速余弦正逆變換(FCT&IFCT)等關鍵代碼，其他代碼如簡單離散傅里葉變換(DFT)與余弦變換(DCT)的代碼請見附件。
1. 主函數與主要參數
constexpr auto AMAX = 100; //調整幅度值可視化參數
constexpr auto CMAX = 50; //調整余弦變換可視化參數
//可選的照片："photo" "Gundam" "lena"
static String name = "Gundam";
//將圖像存儲容器中的值映射到[0,255]的灰度區間內
enum StandardType
{
AMPLITUDE, //幅度圖
PHASE, //相位圖
COSINE, //余弦變換圖
SIFFT, //傅里葉逆變換
SIFCT //余弦逆變換
};
int Shift(Mat& src, Mat& dst); //中心化
int Standard(Container src, Mat& dst, int type); //映射到[0,255]區間
int DFT(Mat src, Container& A, Container& φ); //二維離散傅里葉變換
int DCT(Mat src, Container& C); //二維離散余弦變換
int FFT(Mat src, Container& R, Container& I, Container& A, Container& φ); //快速傅里葉變換
int FCT(Mat src, Container& C); //快速余弦變換
int IFFT(Container R, Container I, Mat& dst); //快速傅里葉逆變換
int IFCT(Container& C, Mat& dst); //快速余弦逆變換
2. 快速傅里葉變換
int FFT(Mat src, Container& R, Container& I, Container& A, Container& φ) {
double M = src.rows; double N = src.cols;
for (int v = 0; v < N; v++) {
vector<double> listR, listI;
for (int x = 0; x < M; x++) {
double r = 0.0, i = 0.0;
for (int y = 0; y < N; y++) {
double t = -2 * PI * v * y / N;
r += src.ptr<uchar>(x)[y] * cos(t);
i += src.ptr<uchar>(x)[y] * sin(t);
}
listR.push_back(r);
listI.push_back(i);
}
for (int u = 0; u < M; u++) {
double r = 0.0, i = 0.0;
for (int x = 0; x < M; x++) {
double t = -2 * PI * u * x / M;
r += listR[x] * cos(t) - listI[x] * sin(t);
i += listR[x] * sin(t) + listI[x] * cos(t);
}
R[u].push_back(r); I[u].push_back(i);
A[u].push_back(sqrt(r * r + i * i));
φ[u].push_back(atan(i / r));
}
}
return 0;
}
3. 快速余弦變換
int FCT(Mat src, Container& C) {
for (int v = 0; v < N; v++) {
vector<double> list;
for (int x = 0; x < M; x++) {
double c = 0.0;
for (int y = 0; y < N; y++) {
double t = (2.0 * y + 1.0) * v * PI / (double)(2.0 * N);
t = (double)src.ptr<uchar>(x)[y] * cos(t);
if (v == 0) t /= sqrt(2);
c += t;
}
list.push_back(c);
}
for (int u = 0; u < M; u++) {
double c = 0.0;
for (int x = 0; x < M; x++) {
double t = (2.0 * x + 1.0) * u * PI / (double)(2.0 * M);
t = list[x] * cos(t);
if (u == 0) t /= sqrt(2);
c += t;
}
C[u].push_back(c * 2 / sqrt(M * N));
}
}
return 0;
}
4. 快速傅里葉逆變換
int IFFT(Container R, Container I, Mat& dst) {
double M = dst.rows; double N = dst.cols; Container Cdst(M);
for (int v = 0; v < N; v++) {
vector<double> listR, listI;
for (int x = 0; x < M; x++) {
double r = 0.0, i = 0.0;
for (int y = 0; y < N; y++) {
double t = 2 * PI * v * y / N;
r += R[x][y] * cos(t) - I[x][y] * sin(t);
i += R[x][y] * sin(t) + I[x][y] * cos(t);
}
listR.push_back(r);
listI.push_back(i);
}
for (int u = 0; u < M; u++) {
double ifft = 0.0;
for (int x = 0; x < M; x++) {
double t = 2 * PI * u * x / M;
ifft += listR[x] * cos(t) - listI[x] * sin(t);
}
Cdst[u].push_back(ifft / (M * N));
}
}
Standard(Cdst, dst, SIFFT);
return 0;
}
5. 快速余弦逆變換
int IFCT(Container& C, Mat& dst) {
double M = dst.rows; double N = dst.cols; Container Cdst(M);
for (int v = 0; v < N; v++) {
vector<double> list;
for (int x = 0; x < M; x++) {
double c = 0.0;
for (int y = 0; y < N; y++) {
double t = (2.0 * y + 1.0) * v * PI / (double)(2.0 * N);
t = C[x][y] * cos(t);
if (v == 0) t /= sqrt(2);
c += t;
}
list.push_back(c);
}
for (int u = 0; u < M; u++) {
double c = 0.0;
for (int x = 0; x < M; x++) {
double t = (2.0 * x + 1.0) * u * PI / (double)(2.0 * M);
t = list[x] * cos(t);
if (u == 0) t /= sqrt(2);
c += t;
}
Cdst[u].push_back(c * 2.0 / sqrt(M * N));
}
}
Standard(Cdst, dst, SIFCT);
return 0;
}
結果展示
1. 基本展示
從左到右依次為 原圖 - 幅度圖 - 余弦變換圖 - 相位圖
2. 傅里葉變換
為了更好的展示兩種變換的特點，避免將數據映射到[0,255]灰度區間所帶來的數據丟失，我在三維空間坐標系中重新繪制了幅度圖與相位圖，如下圖所示。
（1）幅度圖
（2）中心化以后的幅度圖
（3）相位圖
2. 余弦變換
從下圖我們可以看出，余弦變換后圖像信號的能量集中於一角，這是余弦變換最顯著的特點。
作業小結
- 本次作業使用C++語言進行編寫，在編寫過程中，數字精度的控制給我造成了很大的困擾；
- 映射到[0,255]灰度區間時不同映射函數的選擇也給我帶來了很大的麻煩，最后我為不同的圖像設置了不同的映射函數，通過枚舉變量作為參數進行選擇；
- 一直沒有找到合適的可視化方法，直到 Origin 出現在我的面前；
- 傅里葉變換與余弦變換的原理掌握還不夠透徹，需要進一步理論的加強；
- 本想要做濾波等傅里葉變換的實際應用，但終是沒有時間，下次一定，下次一定！