在網上看到很多關於傅立葉變換的內容, 但是沒找到具體工程上完整的一個例子
例如把一個紋理轉化為頻譜圖和相位 然後利用頻譜和相位在轉化回來
於是就自己做一個好了
如果有不對之處請使勁噴
然後如果你比較熟悉只想看工程部分的內容, 可以酌情跳過
先看一眼結果:
放出Git 地址:
https://github.com/lingzerg/LingzergDemo
在看一眼流程:
然後我們慢慢解釋
什麼是信號?先把思路捋清楚
講傅立葉變換的定義有很多, 推薦大家各種搜一下
比如3B1B的傅立葉變換:
https://www.bilibili.com/video/BV1pW411J7s8
比如這篇文章:
https://zhuanlan.zhihu.com/p/19763358
我簡單描述下, 正常情況下 我們看到的信號是這樣的:
這種信號你可以想像是聲音, 在時間軸上連續不斷變換,自然界中的信號我們認為他是連續的, 比如這樣:
橫軸t就是時間, 所以這個叫做時域這就是一個信號在時間軸上不斷變化, 所以這種信號我們叫一個信號的時域表示
然後為了讓計算機可以記錄描述這種信號,我們會這樣:
對信號進行離散採樣以一個固定間隔對信號進行採樣, 得到的數據包含這個信號在當前採樣點的振幅和頻率
公式描述是:
f是頻率(frequency), A是振幅(Amplitude), phase是相位,x是時間, 得到的是振幅
所以我們只要是討論計算機領域的信號, 就一定是離散的
推薦大家去坐標系中手動拉一拉這些變量感受下這些變量的作用:
Desmos | Graphing Calculator
我們拿到一個信號,有時候想過濾掉高頻部分,例如通話降噪,還有一些不想要的低頻部分
這時我們可以用複數坐標系去描述一個信號, 讓信號的另外一個特徵暴露出來,以便於我們處理
複數坐標系:
使用複數平面描述信號有什麼好處呢?
在複數和實數構建的坐標系中,信號的時間周期就變成了一個環, 信號採樣就是不斷圍繞這個圈採樣
模長就是振幅,
那問題就在於,我們如何把一個我們抽樣得到的信號樣本轉換成複數形式?
先把歐拉公式扔這裡:
然後歐拉公式實際上是逆時針旋轉的,而正變換是順時針的
所以我們需要改成
而逆變換就要換回來
接著把頻率 和時間 帶入歐拉公式另一側, 除以 是因為N次採樣
在乘以2π 轉成弧度:
就如同公式裡的描述:
而這個又可以展開寫一篇了, 所以詳細描述我推薦諸位看我上面發的哪兩個連結
最後我們反正知道把一個信號的採樣結果帶入離散傅立葉變換的公式就可以得到一組複數的描述
在這個公式裡 實際上 x 就是時間點, u就是頻率,
小 就是每個單位時間的採樣點採樣得到的振幅,
最後會得到一個複數, 而逆變換就是反過來:
大 就是傅立葉變換得到的複數
二維傅立葉變換?一維傅立葉變換和二維的差別,就在於二維還要在套一層,你需要遍歷整個二維空間
二維的離散逆傅立葉變換:
而在圖像處理中, 二維傅立葉變換上任意一個像素可以認為是一整個平面波:
還可以把圖片上的每個像素看成一個在x軸上方的信號,
比如這樣:
一個sin(x)+1的信號接著將每個採樣點帶入公式就可以得到每個像素對應的複數
我開始做傅立葉變換的時候非常奇怪
我感覺傅立葉變換和逆變換的公式之間並沒有頻域和相位啊
那這個流程到底是什麼?
我先上一張圖然後慢慢解釋
這個基本上就是傅立葉變換的完整流程了
也就是說, 傅立葉變換原始公式裡是只有複數的,但是通過複數我們可以得到每個像素點的膜長, 還有弧度:
你不做這一步, 直接把複數的實部和虛部分別保存在紋理的r和g通道也是完全沒問題的
但是這就失去了對二維圖片進行傅立葉變換的意義
我們只有獲取了頻域和時域,才能根據頻域和時域處理圖片
我看過很多demo, 都是直接保存複數,而沒有對頻域和相位進行獨立的轉換
那我們怎麼在unity裡實現呢?首先我選擇用compute shader來實現這個
因為感覺CS比較帥把
在Cshader中我們創建一個複數的結構體,並寫好常用方法:
struct Complex{ float real; float imaginary;};
Complex getComplex(float num1, float num2){ Complex c; c.real = num1; c.imaginary = num2; return c;}
Complex getAdd(Complex c, Complex o){ c.real = c.real + o.real; c.imaginary = c.imaginary + o.imaginary; return c;}
Complex getSubtract(Complex c,Complex o){ c.real = c.real - o.real; c.imaginary = c.imaginary - o.imaginary; return c;} Complex getMultiply(Complex c,Complex o){ float tmp = c.real * o.real - c.imaginary * o.imaginary; c.imaginary = c.real * o.imaginary + c.imaginary * o.real; c.real = tmp; return c;}
Complex getMultiplyFloat(Complex c,float x){ c.real = c.real *x; c.imaginary = c.imaginary * x; return c;}
Complex getEuler(float x) { Complex c; c.real = cos(x); c.imaginary = sin(x); return c;}
int2 getCoordinate(int x, int y) { x = x - 31; y = y - 31; return int2(x,y);}然後我們修改一下CS的核心, 按照公式做積分,詳細描述在注釋中:
[numthreads(16,16,1)]void Fourier (uint3 id : SV_DispatchThreadID){ float flag = -1 * 2 * 3.14 / TextureSize; int2 coord = id.xy; Complex sumRow; sumRow.real = 0; sumRow.imaginary = 0;
for(int i = 0; i< TextureSize; i++) { Complex complexU = getEuler(flag*i*coord.x); Complex sumColumn; sumColumn.real = 0; sumColumn.imaginary = 0; for(int j = 0; j < TextureSize; j++) { Complex complexV = getEuler(flag*j*coord.y); complexV = getMultiplyFloat(complexV, originalImg[int2(i,j)].r); sumColumn = getAdd(sumColumn,complexV); }
sumRow = getAdd(sumRow,getMultiply(sumColumn, complexU)); } float amplitude = sqrt(sumRow.real * sumRow.real + sumRow.imaginary * sumRow.imaginary); float phase = atan2(sumRow.imaginary , sumRow.real);
rtFourierSpectrum[id.xy] = amplitude/(TextureSize*TextureSize); rtFourierPhase[id.xy] = phase;
rtReal[id.xy] = amplitude * cos(phase); rtImaginary[id.xy] = amplitude * sin(phase); }這時候我們就得到了頻譜圖和相位圖,你可以在這裡對頻譜或者相位進行一些操作
例如我注釋掉那個 //clamp(spectrum/TextureSize,0,0.5); 可以過濾一些高頻信息
或者用各種濾波器在這裡處理頻譜圖, 頻譜圖代表圖片的灰度信息,可以理解為亮度, 而相位代表位置信息.
然後用逆變換還原回去 我們就直接還原回去:
void FourierInverse(uint3 id : SV_DispatchThreadID){ float flag = 2 * 3.14 / TextureSize; int2 coord = id.xy; Complex sumRow; sumRow.real = 0; sumRow.imaginary = 0;
for(int i = 0; i< TextureSize; i++) { Complex complexU = getEuler(flag*i*coord.x); Complex sumColumn; sumColumn.real = 0; sumColumn.imaginary = 0; for(int j = 0; j < TextureSize; j++) { Complex complexV = getEuler(flag*j*coord.y); Complex sampleComplex; float amplitude = rtFourierSpectrum[int2(i,j)].r; float phase = rtFourierPhase[int2(i,j)].r;
sampleComplex.real = amplitude * cos(phase); sampleComplex.imaginary = amplitude * sin(phase);; complexV = getMultiply(complexV,sampleComplex); sumColumn = getAdd(sumColumn,complexV); } sumRow = getAdd(sumRow,getMultiply(sumColumn, complexU)); } float value = sqrt(sumRow.real * sumRow.real + sumRow.imaginary * sumRow.imaginary); rtFourierInverse[id.xy] = value; }最後給一點tips:
聲明：發布此文是出於傳遞更多知識以供交流學習之目的。若有來源標註錯誤或侵犯了您的合法權益，請作者持權屬證明與我們聯繫，我們將及時更正、刪除，謝謝。
作者：lingzerg
原文：https://zhuanlan.zhihu.com/p/192659303
More：【微信公眾號】 u3dnotes