快速傅里叶变换 java_二维快速傅里叶变换的java实现
圖像處理與模式識別課作業(yè)，沒學(xué)過信號與系統(tǒng)(哭暈)。
惡補(bǔ)了兩天岡薩里斯的書，寫一下實(shí)現(xiàn)原理及過程
看了網(wǎng)絡(luò)上很多版本的概念講解，道理我都懂，但是在將算法遷移到傅里葉變換的實(shí)現(xiàn)上時(shí)，出現(xiàn)了一些問題。接下來我會簡要介紹快速傅里葉變換是如何加速的，著重寫如何將加速算法應(yīng)用到傅里葉變換上。
二維快速傅里葉變換原理介紹
1.1普通的二維傅里葉變換
?二維傅里葉變換的公式如下：
\[F(u, v) = \sum^{M-1}_{x=0}\sum^{N-1}_{y=0}f(x, y)e^{-j2\pi({ux}/M+{vy}/N)}
\]
對上述表達(dá)式調(diào)整和分順序，先對x求和，再對y求和。我們可以將二維離散傅里葉變換公式轉(zhuǎn)化為兩次一維離散傅里葉變換。一維離散傅里葉變換公式如下：
\[F(u)=\sum^{N-1}_{x=0}f(x)e^{{{-j2\pi}ux}/N}\quad u = 0,1,2...N-1
\]
?從變換公式中我們可以看出來，由于我們要求N個(gè)\(F(u)\)，每個(gè)\(F(u)\)需要進(jìn)行N次運(yùn)算。使用這樣算法的時(shí)間復(fù)雜度是\(O(n^2)\)的。這樣的算法在N特別大時(shí)，會比較慢。在這里引入快速傅里葉變換的加速算法。
1.2快速傅里葉變換原理
?相信大家在看過別的博客內(nèi)容后已經(jīng)掌握了快速傅里葉變換的原理。我在這里就簡要再提一下。快速傅里葉變換其中一個(gè)重要的依據(jù)就是函數(shù)的點(diǎn)值對表示。其實(shí)我們對離散函數(shù)求得的傅里葉變換的結(jié)果就是一個(gè)點(diǎn)值對表示(畢竟我們沒有求出傅里葉變換函數(shù)的系數(shù))。
?一個(gè)n次多項(xiàng)式函數(shù)，注意：這里的n是2的整數(shù)冪，如果一個(gè)函數(shù)的最高次沒有達(dá)到2的整數(shù)冪，我們可以認(rèn)為此項(xiàng)系數(shù)為0，不影響函數(shù)的值
\[A(x) = a_0 + a_1x+a_2x^2+...+a_{n-1}x^{n-1}
\]
在奇偶項(xiàng)合并后，我們可以得到
\[A(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x+a_3x^3+...+a_{n-1}x^{n-1})\\
A(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})\\
A(x)=A_{even}(x^2)+xA_{odd}(x^2)
\]
?我們將\(w^k_n\quad k\in(0, 1, 2,...,n-1)\)帶入A(x)。可以得到A(x)的點(diǎn)值對表示。但是這種點(diǎn)值對的計(jì)算方法，仍然是\(O(n^2)\)的。如何將計(jì)算加速呢？這里就用到我們?nèi)∏傻腬(w^k_n\)上了。
?簡要介紹一下\(w^k_n\)，它是正整數(shù)1的n次單位根。具體值為：
\[w^k_n=e^{i{2\pi}k/n}=cos(2\pi k/n)+i*sin(2\pi k/n)
\]
從虛數(shù)相乘等于在圓上的旋轉(zhuǎn)，我們可以看出它確實(shí)是1的n次單位根的集合。它有幾個(gè)特殊的性質(zhì)
\[w^{k+n/2}_n=e^{i2\pi (k+n/2)/n}=e^{i2\pi k/n}*e^{i\pi}=e^{i2\pi k/n}*(cos(\pi)+i*sin(\pi))=-e^{i2\pi k/n}=-w^k_n\\
w^{k+n}_n=e^{i2\pi k/n}*e^{i2\pi}=e^{i2\pi k/n}=w^k_n
\]
所以我們在計(jì)算點(diǎn)值對的時(shí)候，有這樣一種關(guān)系：
\[A(w^k_n)=A_{even}((w^k_n)^2)+w^k_n*A_{odd}((w^k_n)^2)
-->A(w^k_n)=A_{even}(w^{2k}_n)+w^k_n*A_{odd}(w^{2k}_n)\\
A(w^{k+n/2}_n)=A_{even}((w^{k+n/2}_n)^2)+w^{k+n/2}_n*A_{odd}((w^{k+n/2}_n)^2)-->
A(w^{k+n/2}_n)=A_{even}(w^{2k+n}_n)-w^k_n*A_{odd}(w^{2k+n}_n)-->A(w^k_n)=A_{even}(w^{2k}_n)-w^k_n*A_{odd}(w^{2k}_n)
\]
整理一下，結(jié)果如下
\[A(w^k_n)=A_{even}(w^{2k}_n)+w^k_n*A_{odd}(w^{2k}_n)\\
A(w^{k+n/2}_n)=A_{even}(w^{2k}_n)-w^k_n*A_{odd}(w^{2k}_n)
\]
則我們求出\(A_{even}(w^{2k}_n)\)和\(A_{odd}(w^{2k}_n)\)就可以得到兩個(gè)值。根據(jù)我們算法課學(xué)的分治+遞歸的思想。遞歸層數(shù)為logn次，每層計(jì)算時(shí)間復(fù)雜度為\(O(n)\)，則總的時(shí)間復(fù)雜度為\(O(nlogn)\)。這樣我們就做到了快速計(jì)算。
1.3將算法應(yīng)用到傅里葉變換
?之前我們大致了解了傅里葉變換的原理，那么如何將這個(gè)算法加速我們傅里葉變換的計(jì)算呢？一維離散傅里葉變換的公式為：
\[F(u)=\sum^{N-1}_{x=0}f(x)e^{{{-j2\pi}ux}/N}\quad u = 0,1,2...N-1
\]
是不是看到了\(e^{{{-j2\pi}ux}/N}\)? 注意：我們接下來使用的\(w^k_n\)=\(e^{-j2\pi k/n}\)。可以證明，這并不影響我們之前介紹的性質(zhì)。
?則我們的傅里葉變換公式可以寫成
\[F(u)=\sum^{N-1}_{x=0}f(x)w^{ux}_n
\]
那么我們應(yīng)用加速的地方就出現(xiàn)了
\[F(u)=\sum_{t=0}^{\frac{n}{2}-1}f_{even}(t)w^{2ut}_n+w^u_n * \sum_{t=0}^{\frac{n}{2}-1}f_{odd}(t)w^{2ut}_n\\
F(u+n/2)=\sum_{t=0}^{\frac{n}{2}-1}f_{even}(t)w^{2ut}_n-w^u_n * \sum_{t=0}^{\frac{n}{2}-1}f_{odd}(t)w^{2ut}_n
\]
是不是和我們之前看到的算法原理的地方很像？是的，我們快速傅里葉變換就是利用這個(gè)公式做的。利用這個(gè)公式，一次就可以計(jì)算出兩個(gè)值，加速了傅里葉變換計(jì)算的層數(shù)。同樣的，這個(gè)算法共有l(wèi)ogn層，每層的時(shí)間復(fù)雜度為\(O(n)\)。整體算法的時(shí)間復(fù)雜度為\(O(nlogn)\)。
這里還有一個(gè)trick，即蝶形算法。我們在這里舉個(gè)蝶形算法的例子。假設(shè)有下列一組數(shù)
\[f(0)\quad f(1)\quad f(2)\quad f(3)\quad f(4)\quad f(5)\quad f(6)\quad f(7)
\]
第一次劃分奇偶后，我們得到
\[[f(0)\quad f(2)\quad f(4)\quad f(6)]\quad [f(1)\quad f(3)\quad f(5)\quad f(7)]
\]
第二次劃分奇偶后，結(jié)果
\[[f(0)\quad f(4)]\quad [f(2)\quad f(6)]\quad [f(1)\quad f(5)]\quad [f(3)\quad f(7)]
\]
第三次到遞歸終點(diǎn)，我們看一下上訴編號的二進(jìn)制碼
\[000\quad 100\quad 010\quad 110\quad 001\quad 101\quad 011\quad 111
\]
將二進(jìn)制反序
\[000\quad 001\quad 010\quad 011\quad 100\quad 101\quad 110\quad 111
\]
是不是就成了一個(gè)自然數(shù)列？我們可以根據(jù)這個(gè)蝶形的方式，先將\(f(x)\)排列好，然后以循環(huán)的方式計(jì)算最終的結(jié)果。
核心java代碼
import java.util.Comparator;
import java.util.ArrayList;
public class FFT {
private ArrayList data;
private ArrayList buff;
public FFT(ArrayList list, int pow){
data = new ArrayList<>();
for(int i = 0; i < list.size(); i++){
sortData tmp = new sortData(list.get(i), reverse(i, pow));
data.add(tmp);
}
Comparator comparator = (o1, o2) -> Integer.compare(o2.getRevNum(), o1.getRevNum());
data.sort(comparator);
buff = new ArrayList<>();
}
public ArrayList trans(){
recursion(data.size());
return buff;
}
private void recursion(int half){
for (sortData tmpData : data) {
// init the array
Complex tmp = new Complex(1.0);
tmp.multiply(tmpData.getValue());
buff.add(tmp);
}
int count = 2;
while(half >= count) {
for (int i = 0; i < data.size(); i = i + count) {
for (int j = i; j < i + count / 2; j++) {
Complex tmp1 = buff.get(j);
Complex tmp2 = buff.get(j + count / 2);
Complex tmpComplex = tmp1.complexClone();
Complex w = root(count, 1, j - i);
tmp2.multiply(w);
tmpComplex.add(tmp2);
tmp1.subtract(tmp2);
buff.set(j, tmpComplex);
buff.set(j + count / 2, tmp1);
}
}
count = count * 2;
}
}
private int reverse(int value, int pow){
int res = Integer.reverse(value);
res = res >> (Integer.SIZE - pow);
res = res & ((1 << pow) - 1);
return res;
}
private Complex root(int k, int x, int u){
double real = Math.cos(Math.PI * 2 * u * x / k);
double imaginary = -1 * Math.sin(Math.PI * 2 * u * x / k);
return new Complex(real, imaginary);
}
}
class sortData{
private Complex value;
private int revNum;
public sortData(Complex o, int revNum){
value = o;
this.revNum = revNum;
}
public int getRevNum(){
return revNum;
}
public Complex getValue(){
return value;
}
}
recursion函數(shù)是計(jì)算構(gòu)造方法中傳入的ArrayList的傅里葉變換，計(jì)算結(jié)果存在buff中。
reverse函數(shù)是利用蝶形算法的原理，計(jì)算出每個(gè)采樣點(diǎn)二進(jìn)制的反轉(zhuǎn)結(jié)果。
root函數(shù)是計(jì)算1的n次單位結(jié)果\(w^{ux}_k\)。
Complex是一個(gè)復(fù)數(shù)類，具體內(nèi)容如下。
public class Complex {
private double real;
private double imaginary;
public Complex(){
real = 0;
imaginary = 0;
}
public Complex(double num){
real = num;
imaginary = 0;
}
public Complex(double num1, double num2){
real = num1;
imaginary = num2;
}
public Complex complexClone(){
return new Complex(real, imaginary);
}
public void add(Complex o){
real = real + o.real;
imaginary = imaginary + o.imaginary;
}
public void subtract(Complex o){
real = real - o.real;
imaginary = imaginary - o.imaginary;
}
public void multiply(Complex o){
double tmp = real * o.real - imaginary * o.imaginary;
imaginary = real * o.imaginary + imaginary * o.real;
real = tmp;
}
public void absLog(){
real = Math.sqrt(real * real + imaginary * imaginary)/50;
if(real > 255){
real = 255;
}
imaginary = 0;
}
public int getReal(){
return (int)real;
}
@Override
public String toString(){
if (imaginary == 0) return real + "";
if (real == 0) return imaginary + "i";
if (imaginary < 0) return real + " - " + (-imaginary) + "i";
return real + " + " + imaginary + "i";
}
}
填坑完畢，轉(zhuǎn)發(fā)需注明出處。
總結(jié)
以上是生活随笔為你收集整理的快速傅里叶变换 java_二维快速傅里叶变换的java实现的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： Bear and Big Brother
- 下一篇： 关于模拟题的一些弱鸡总结