您的位置:首页 > 编程语言 > C语言/C++

最近开始学习图像识别,第一步FFT

2007-09-04 10:14 603 查看
代码根据csdn论坛上一位达人的C++代码改编:

public class FFT {

/** Creates a new instance of FFT */
public FFT(int _ex) {

ex = _ex;
N = (int) Math.pow(2, ex);
omegaRe = new double
;
omegaIm = new double
;
for (int k = 0; k < N; k++) {
omegaRe[k] = Math.cos(2 * Math.PI / N * k);
omegaIm[k] = Math.sin(-2 * Math.PI / N * k);
}
}

double omegaRe[];
double omegaIm[];
int ex;
int N;

public int rev(int x) {
int in = x;
int ret = 0;
for (int i = 0; i < ex; i++) {
ret = ret | (in % 2 << ex - i - 1);
in = in >> 1;
}
return ret;
}

public void fft(double[] inputRe, double[] inputIm, double[] outputRe,
double[] outputIm) {

assert (inputRe.length >= N);
assert (inputRe.length == inputIm.length);
assert (outputRe.length == outputIm.length);
assert (inputRe.length == outputRe.length);

int k; // 当前树的深度
int j;

int nBtFlyLen = 0;

if (inputIm == null) {
inputIm = new double[inputRe.length];
}

// 变换需要的工作空间
double[] work1Re = new double
;
double[] work1Im = new double
;
double[] work2Re = new double
;
double[] work2Im = new double
;

// 临时变量
double[] tmpRe;
double[] tmpIm;
// 初始化,写入数据
System.arraycopy(inputRe, 0, work1Re, 0, N);
System.arraycopy(inputIm, 0, work1Im, 0, N);
// 临时变量
int nInter = 0;

// 蝶形算法进行快速傅立叶变换
for (k = 0; k < ex; k++) {
for (j = 0; j < (1 << k); j++) {
// 计算长度
nBtFlyLen = 1 << (ex - k);
// 倒序重排,加权计算
for (int i = 0; i < nBtFlyLen / 2; i++) {
nInter = j * nBtFlyLen;
work2Re[i + nInter] = work1Re[i + nInter]
+ work1Re[i + nInter + nBtFlyLen / 2];
work2Im[i + nInter] = work1Im[i + nInter]
+ work1Im[i + nInter + nBtFlyLen / 2];

double a = work1Re[i + nInter]
- work1Re[i + nInter + nBtFlyLen / 2];
double b = work1Im[i + nInter]
- work1Im[i + nInter + nBtFlyLen / 2];
double c = omegaRe[i * (1 << k)];
double d = omegaIm[i * (1 << k)];
work2Re[i + nInter + nBtFlyLen / 2] = a * c - b * d;
work2Im[i + nInter + nBtFlyLen / 2] = a * d + b * c;
}
}
// 交换 work1和work2的数据
tmpRe = work1Re;
tmpIm = work1Im;
work1Re = work2Re;
work1Im = work2Im;
work2Re = tmpRe;
work2Im = tmpIm;
}

// 重新排序
for (j = 0; j < N; j++) {
nInter = 0;
for (int i = 0; i < ex; i++) {
if ((j & (1 << i)) != 0) {
nInter += 1 << (ex - i - 1);
}
}
outputRe[j] = work1Re[nInter];
outputIm[j] = work1Im[nInter];
}

}

public void ifft(double[] inputRe, double[] inputIm, double[] outputRe,
double[] outputIm) {

assert (inputRe.length >= N);
assert (inputRe.length == inputIm.length);
assert (outputRe.length == outputIm.length);
assert (inputRe.length == outputRe.length);

double[] transformRe = new double
;
double[] transformIm = new double
;

for (int i = 0; i < N; i++) {
transformRe[i] = inputRe[i];
transformIm[i] = -inputIm[i];
}
fft(transformRe, transformIm, outputRe, outputIm);
for (int i = 0; i < N; i++) {
outputRe[i] = outputRe[i] / N;
outputIm[i] = -outputIm[i] / N;
}
}

public void fft_2d(double[][] inputRe, double[][] inputIm,
double[][] outputRe, double[][] outputIm) {

// 检查不强, 没有循环检查二维数组的第二维
assert (inputRe.length >= N);
assert (inputRe.length == inputIm.length);
assert (outputRe.length == outputIm.length);
assert (inputRe.length == outputRe.length);

if (inputIm == null)
inputIm = new double

;

for (int i = 0; i < N; i++) {
fft(inputRe[i], inputIm[i], outputRe[i], outputIm[i]);
}

// 重新为in分配一块空间
double[][] inRe = new double

;
double[][] inIm = new double

;

// 转置
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
inRe[i][j] = outputRe[j][i];
inIm[i][j] = outputIm[j][i];
}
}
for (int i = 0; i < N; i++) {
fft(inRe[i], inIm[i], outputRe[i], outputIm[i]);
}
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
inRe[i][j] = outputRe[j][i];
inIm[i][j] = outputIm[j][i];
}
}
for (int i = 0; i < N; i++) {
System.arraycopy(inRe[i], 0, outputRe[i], 0, N);
System.arraycopy(inIm[i], 0, outputIm[i], 0, N);
}

}

public void ifft_2d(double[][] inputRe, double[][] inputIm,
double[][] outputRe, double[][] outputIm) {

double[][] workRe = new double

;
double[][] workIm = new double

;

for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
workRe[i][j] = inputRe[i][j];
workIm[i][j] = -inputIm[i][j];
}
}
fft_2d(workRe, workIm, outputRe, outputIm);
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
outputRe[i][j] = outputRe[i][j] / (N * N);
outputIm[i][j] = -outputIm[i][j] / (N * N);
}
}
}
}

[/code]

创建FFT对象的时候指定层数ex,也就是点数N= pow(2,ex); 构造函数中会计算omega(0到N-1)节约时间。
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签:  fft im null c++ class 算法