頻域處理：傅立葉變換及小波變換
頻域處理：傅立葉變換及小波變換
引言
影象處理–>頻域處理–>傅立葉變換、小波變換。用另一種方法來觀察世界的話，你會發現世界是永恆不變的。
1、傅立葉變換
尤拉公式所描繪的，是一個隨著時間變化，在複平面上做圓周運動的點，隨著時間的改變，在時間軸上就成了一條螺旋線。如果只看它的實數部分，也就是螺旋線在左側的投影，就是一個最基礎的餘弦函式。而右側的投影則是一個正弦函式。
由尤拉公式知道：正弦波的疊加
複數就是一個能把平面進行均勻縮放和旋轉的乘子。
通過時域到頻域的變換，得到了一個從側面看的頻譜，但是這個頻譜並沒有包含時域中全部的資訊。因為頻譜只代表每一個對應的正弦波的振幅是多少，而沒有提到相位。基礎的正弦波A.sin(wt+θ)中，振幅，頻率，相位缺一不可，不同相位決定了波的位置，所以對於頻域分析，僅僅有頻譜（振幅譜）是不夠的，還需要一個相位譜。
糾正一個概念：時間差並不是相位差。如果將全部週期看作2Pi或者360度的話，相位差則是時間差在一個週期中所佔的比例。我們將時間差除週期再乘2Pi，就得到了相位差。
在完整的立體圖中，將投影得到的時間差依次除以所在頻率的週期，就得到了最下面的相位譜。所以，頻譜是從側面看，相位譜是從下面看。
短時傅立葉變換（Short-time Fourier Transform,STFT）也叫加窗傅立葉變換，顧名思義，就是因為傅立葉變換的時域太長了，所以要弄短一點，這樣就有了區域性性。
定義：把整個時域過程分解成無數個等長的小過程，每個小過程近似平穩，再傅立葉變換，就知道在哪個時間點上出現了什麼頻率了。”這就是短時傅立葉變換。
時域上分成一段一段做FFT，就可以知道頻率成分隨著時間的變化情況。
2、小波變換
小波變換
對於加窗傅立葉變換讓人頭疼的就是視窗的大小問題，如果讓視窗的大小可以改變，就完美了。小波就是基於這個思路，但是不同的是。STFT是給訊號加窗，分段做FFT；而小波變換並沒有採用窗的思想，更沒有做傅立葉變換。小波直接把傅立葉變換的基給換了——將無限長的三角函式基換成了有限長的會衰減的小波基。這樣不僅能夠獲取頻率，還可以定位到時間。即是做傅立葉變換隻能得到一個頻譜，做小波變換卻可以得到一個時頻譜。
從公式可以看出，不同於傅立葉變換，變數只有頻率ω，小波變換有兩個變數：尺度a（scale）和平移量 τ（translation）。尺度a控制小波函式的伸縮，平移量 τ控制小波函式的平移。尺度就對應於頻率（反比），平移量 τ就對應於時間。
當伸縮、平移到這麼一種重合情況時，也會相乘得到一個大的值。這時候和傅立葉變換不同的是，這不僅可以知道訊號有這樣頻率的成分，而且知道它在時域上存在的具體位置。而當每個尺度下都平移著和訊號乘過一遍後，就知道訊號在每個位置都包含哪些頻率成分，就可以對非穩定訊號做時頻分析。
3、程式
fft:
import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt
img = cv.imread('./opencv.png',0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20*np.log(np.abs(fshift))
plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()
#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>
using namespace std;
using namespace cv;
int main()
{
Mat I = imread("lena.jpg", IMREAD_GRAYSCALE); //讀入影象灰度圖
//判斷影象是否載入成功
if (I.empty())
{
cout << "影象載入失敗!" << endl;
return -1;
}
else
cout << "影象載入成功!" << endl << endl;
Mat padded; //以0填充輸入影象矩陣
int m = getOptimalDFTSize(I.rows);
int n = getOptimalDFTSize(I.cols);
//填充輸入影象I，輸入矩陣為padded，上方和左方不做填充處理
copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));
Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F) };
Mat complexI;
merge(planes, 2, complexI); //將planes融合合併成一個多通道陣列complexI
dft(complexI, complexI); //進行傅立葉變換
//計算幅值，轉換到對數尺度(logarithmic scale)
//=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
split(complexI, planes); //planes[0] = Re(DFT(I),planes[1] = Im(DFT(I))
//即planes[0]為實部,planes[1]為虛部
magnitude(planes[0], planes[1], planes[0]); //planes[0] = magnitude
Mat magI = planes[0];
magI += Scalar::all(1);
log(magI, magI); //轉換到對數尺度(logarithmic scale)
//如果有奇數行或列，則對頻譜進行裁剪
magI = magI(Rect(0, 0, magI.cols&-2, magI.rows&-2));
//重新排列傅立葉影象中的象限，使得原點位於影象中心
int cx = magI.cols / 2;
int cy = magI.rows / 2;
Mat q0(magI, Rect(0, 0, cx, cy)); //左上角影象劃定ROI區域
Mat q1(magI, Rect(cx, 0, cx, cy)); //右上角影象
Mat q2(magI, Rect(0, cy, cx, cy)); //左下角影象
Mat q3(magI, Rect(cx, cy, cx, cy)); //右下角影象
//變換左上角和右下角象限
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
//變換右上角和左下角象限
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
//歸一化處理，用0-1之間的浮點數將矩陣變換為可視的影象格式
normalize(magI, magI, 0, 1, CV_MINMAX);
imshow("輸入影象", I);
imshow("頻譜圖", magI);
waitKey(0);
return 0;
}