FFT任意点变换源码(1维,2维图像FFT,图像频域相关,频谱相乘)

#include "fft.h" #include <stdio.h> #include <string.h> #include <math.h> #include <stdlib.h> static int pow_int(int x, int y) { return pow(x, y); int m = 1; for (int i = 0; i < y; i++) m *= x; return m; } void spk_fft::FFT1D(spk_complex *src, int n) { spk_complex *buf1 = new spk_complex[n]; spk_complex *buf2 = new spk_complex[n]; memcpy(buf1, src, sizeof(spk_complex)*n); //将arrayBuf数组元素基2抽取并重新排列 //若0、1、2、3、4、5、6、7八点序列对调后变作0、4、2、6、1、5、3、7 for (int r = 1; pow_int(2, r) < n; r++) { int t1 = pow_int(2, r); int t2 = pow_int(2, r - 1); for (int k = 0; k < t1; k++) { for (int i = 0; i < n / t1; i++) { buf2[k*n / t1 + i].real = buf1[k / 2 * n / t2 + i * 2 + k % 2].real; buf2[k*n / t1 + i].imag = buf1[k / 2 * n / t2 + i * 2 + k % 2].imag; } } memcpy(buf1, buf2, sizeof(spk_complex)*n); } //采用蝶型算法进行快速傅立叶变换 //buf1是第r级的输入，buf2存放第r级的输出 for (int r = 1; pow_int(2, r) <= n; r++) { int t1 = pow_int(2, r); for (int k = 0; k < n / t1; k++) { for (int i = t1 / 2; i < t1; i++) { cType c = cos(-2 * PI*(i - t1 / 2) / t1);//加权因子 cType s = sin(-2 * PI*(i - t1 / 2) / t1); buf1[k*t1 + i].real = buf2[k*t1 + i].real*c - buf2[k*t1 + i].imag*s; buf1[k*t1 + i].imag = buf2[k*t1 + i].imag*c + buf2[k*t1 + i].real*s; } } for (int k = 0; k < n / t1; k++) { for (int i = 0; i < t1 / 2; i++) { buf2[k*t1 + i].real = buf1[k*t1 + i].real + buf1[k*t1 + i + t1 / 2].real; buf2[k*t1 + i].imag = buf1[k*t1 + i].imag + buf1[k*t1 + i + t1 / 2].imag; } for (int i = t1 / 2; i < t1; i++) { buf2[k*t1 + i].real = buf1[k*t1 + i - t1 / 2].real - buf1[k*t1 + i].real; buf2[k*t1 + i].imag = buf1[k*t1 + i - t1 / 2].imag - buf1[k*t1 + i].imag; } } memcpy(buf1, buf2, sizeof(spk_complex)*n);//第r级的输出存入buf1,作为下一级的输入数据 } memcpy(src, buf2, sizeof(spk_complex)*n); delete []buf2; delete []buf1; } void spk_fft::IFFT1D(spk_complex *src, int n) { for (int i = 0; i < n; i++) src[i].imag = -src[i].imag; FFT1D(src, n); for (int i = 0; i < n; i++) { src[i].real = src[i].real / n; src[i].imag = -src[i].imag / n; } } void spk_fft::FFT2D(spk_complex* src, int width, int height, spk_complex *&dst) { spk_complex *array = new spk_complex[height]; if (dst == NULL) { dst = new spk_complex[height*width]; } //先纵向一维快速傅立叶变换 for (int u = 0; u < width; u++) { for (int v = 0; v < height; v++) { array[v].real = src[v*width + u].real; array[v].imag = src[v*width + u].imag; } FFT1D(array, height); for (int v = 0; v < height; v++) { dst[v*width+u].real = array[v].real; dst[v*width+u].imag = array[v].imag; } } delete []array; //再横向一维快速傅立叶变换 for (int v = 0; v < height; v++) { FFT1D(dst + v * width, width); } } void spk_fft::IFFT2D(spk_complex * src, int width, int height, spk_complex *& dst) { //先纵向傅立叶反变换 spk_complex *array = new spk_complex[height]; for (int u = 0; u < width; u++) { for (int v = 0; v < height; v++) { array[v].real = src[v*width + u].real; array[v].imag = src[v*width + u].imag; } IFFT1D(array, height); for (int v = 0; v < height; v++) { src[v*width + u].real = array[v].real; src[v*width + u].imag = array[v].imag; } } delete []array; //再横向傅立叶反变换 for (int v = 0; v < height; v++) { IFFT1D(src + v * width, width); } if (dst == NULL) { dst = new spk_complex[width*height]; } for (int i = 0; i < width*height; i++) { dst[i].imag = src[i].imag; dst[i].real = src[i].real; } } int spk_fft::get_size(int x) { int base = 1; for (int i = 0; i < 32; i++) { base = 1 << i; if (x <= base) return base; } return 0; } spk_fft::spk_fft(void) { } /*频域相乘*/ spk_complex* spk_fft::mul(spk_complex *srcA, spk_complex *srcB, int len) { spk_complex *dst = new spk_complex[len]; for (int i = 0; i < len; i++) { dst[i].imag = srcA[i].imag * srcB[i].real - srcA[i].real * srcB[i].imag; dst[i].real = srcA[i].real * srcB[i].real + srcA[i].imag * srcB[i].imag; } return dst; } void spk_fft::fft(cType *src, int _src_row, int _src_col, spk_complex *&dst, int _dst_row, int _dst_col) { int m = _dst_row; int n = _dst_col; int cnt = m * n; spk_complex *src_pad = new spk_complex[cnt]; memset(src_pad, 0, sizeof(spk_complex)*cnt); dst = new spk_complex[cnt]; for (int i = 0; i < _src_row; i++) for (int j = 0; j < _src_col; j++) { src_pad[i*n + j].real = src[i*_src_col+j]; } FFT2D(src_pad, n, m, dst); delete[] src_pad; } void spk_fft::ifft(spk_complex *src, int _src_row, int _src_col, spk_complex *&dst, int _dst_row, int _dst_col) { spk_complex *ifft_dst = NULL; IFFT2D(src, _src_col, _src_row, ifft_dst); if (dst == NULL) { dst = new spk_complex[_dst_row*_dst_col]; } for (int i = 0; i < _dst_row; i++) for (int j = 0; j < _dst_col; j++) { dst[i*_dst_col+j].imag = ifft_dst[i*_src_col+j].imag; dst[i*_dst_col+j].real = ifft_dst[i*_src_col+j].real; } delete[] ifft_dst; } cType * spk_fft::fft_relate( cType *srcA, int Arow, int Acol, cType *srcB, int Brow, int Bcol ) { int m = max(Arow, Brow); int n = max(Acol, Bcol); int _m = get_size(m); int _n = get_size(n); spk_complex *dst = NULL; spk_complex *fftA = NULL; spk_complex *fftB = NULL; spk_complex *mulAB = NULL; fft(srcA, Arow, Acol, fftA, _m,_n); /*对图像A进行FFT变换*/ fft(srcB, Brow, Bcol, fftB, _m,_n); /*对图像B进行FFT变换*/ mulAB = mul(fftA, fftB, _m*_n); /*频域相乘求相关*/ m = m - min(Arow, Brow) + 1; n = n - min(Acol, Bcol) + 1; ifft(mulAB,_m,_n,dst,m,n); /*反变换到时域*/ cType *_dst = new cType[m*n]; for (int i = 0; i < m*n; i++) { _dst[i] = dst[i].real; /*取实部作为相关结果*/ } delete []fftA; delete []fftB; delete []mulAB; delete []dst; return _dst; } /* void fft_test(Mat img) { spk_fft fft; spk_complex *dst_fft = NULL; Mat dstImg; dst_fft = fft.fft(img); dstImg = fft.ifft(dst_fft); for (int j = 0; j < img.rows*img.cols; j++) { if (dstImg.data[j] != img.data[j]) { printf("Err: src[%d]=%d,dst[%d]=%d\r\n", j, img.data[j], j, dstImg.data[j]); } } imshow("原始图像", img); imshow("FFT图像", dstImg); waitKey(0); delete[]dst_fft; } spk_complex * spk_fft::fft(Mat img) { int m = get_size(img.rows); int n = get_size(img.cols); int cnt = m * n; spk_complex *src = new spk_complex[cnt]; memset(src, 0, sizeof(spk_complex)*cnt); spk_complex *dst = new spk_complex[cnt]; //空余部分填充0 for (int i = 0; i < img.rows; i++) for (int j = 0; j < img.cols; j++) { int src_index = i * img.cols + j; int dst_index = i * n + j; src[dst_index].real = *(img.ptr<uchar>(i) + j); } FFT2D(src, n, m, dst); delete []src; src_col = img.cols; src_row = img.rows; dst_row = m; dst_col = n; return dst; } Mat spk_fft::ifft(spk_complex *src, Mat *map) { spk_complex *ifft_dst = NULL; Mat dst = Mat::zeros(src_row,src_col,CV_8UC1); IFFT2D(src, dst_col, dst_row, ifft_dst); for (int i = 0; i < dst.rows; i++) for (int j = 0; j < dst.cols; j++) { int src_index = i*dst_col+j; *(dst.ptr<uchar>(i) + j) = ifft_dst[src_index].real + 0.5; } delete[] ifft_dst; //输出频谱图 if (map != NULL) { Mat fft_map = Mat::zeros(dst_row, dst_col, CV_8UC1); cType t; int i0, j0; for (int i = 0; i < dst_row; i++) for (int j = 0; j < dst_col; j++) { //将频谱移动到图像中央 if (i < dst_row / 2) i0 = i + dst_row / 2; else i0 = i - dst_row / 2; if (j < dst_col / 2) j0 = j + dst_col / 2; else j0 = j - dst_col / 2; int index = i0 * dst_col + j0; t = src[index].real*src[index].real; t += src[index].imag*src