C++实现二维fft和ifft

我是搜遍了没找到实现方法，要么是opencv，要么是dft，然而dft效率差的惊人，随便一个图像都跑不出来

二维傅里叶变换相当于先按行变换，再按列变换

所以首先是一维fft的封装，但封装fft前，得先是复数类的封装

complex.h

#ifndef COMPLEX_H
#define COMPLEX_H
#define MAX_MATRIX_SIZE 4194304 // 2048 * 2048
#define PI 3.141592653
#define EXP 2.718281828
#include
#include
#include
class Complex
{
private:
public:
double real;
double imag;
Complex();
Complex(double rl, double im);
~Complex(void);
// 重载四则运算符号
inline Complex operator +(const Complex &c) {
return Complex(real + c.real, imag + c.imag);
}
inline Complex operator -(const Complex &c) {
return Complex(real - c.real, imag - c.imag);
}
inline Complex operator *(const Complex &c) {
return Complex(real*c.real - imag*c.imag, imag*c.real + real*c.imag);
}
inline Complex operator /(const Complex &c) {
if ((0==c.real) && (0==c.imag)) {
qDebug()<<"11111 ComplexNumber ERROR: divider is 0!";
return Complex(real, imag);
}
return Complex((real*c.real + imag*c.imag) / (c.real*c.real + c.imag*c.imag),
(imag*c.real - real*c.imag) / (c.real*c.real + c.imag*c.imag));
}
inline Complex operator /(const double &c) {
if (0==c) {
qDebug()<<"11111 ComplexNumber ERROR: divider is 0!";
return Complex(real, imag);
}
return Complex(real/c,imag/c);
}
inline Complex operator =(const Complex &c){
real=c.real;
imag=c.imag;
return *this;
}
inline bool operator ==(const Complex &c){
return (real==c.real)&&(imag==c.imag);
}
friend QDebug operator <<(QDebug debug, const Complex &c) {
debug <<"("<<((abs(c.real)<0.000001)?0:c.real) <<"+" <<((abs(c.imag)<0.000001)?0:c.imag)<<"i)";
return debug;
}
double mod(){
return sqrt(real*real+imag*imag);
}
void SetValue(double rl, double im);
};
#endif // COMPLEX_H


complex.cpp

#include "complex.h"
Complex::Complex()
{
real = 0;
imag = 0;
}
Complex::Complex(double rl, double im)
{
real = rl;
imag = im;
}
Complex::~Complex(void)
{
}
void Complex::SetValue(double rl, double im) {
real = rl;
imag = im;
}


然后是一维fft和ifft

dft_one.h

#ifndef DFT_ONE_H
#define DFT_ONE_H
#include
#include "complex.h"
class DFT_one
{
public:
DFT_one(void);
~DFT_one(void);
bool fft(double vec[], int len); // 一维离散傅里叶变换
bool fft(Complex vec[], int len); // 一维离散傅里叶变换
bool ifft(double vec[], int len); // 一维离散傅里叶变换
bool ifft(Complex vec[], int len); // 一维离散傅里叶变换
double *idft(int &ilen); // 一维离散傅里叶逆变换
bool has_dft_vector(); // 是否已存有变换结果
void clear_dft_vector(); // 清除已有的变换结果
void clear_idft_vector();
void print(); // 打印变换结果
Complex *m_dft_vector; // 保存变换结果的容器
Complex *m_idft_vector;
QVectord1_vector;
QVectorid1_vector;
private:
bool m_has_dft_vector;
int m_dft_vector_size; // 变换结果的长度
bool m_has_idft_vector;
int m_idft_vector_size;
};
#endif // DFT_ONE_H


dft_one.cpp

#include "dft_one.h"
DFT_one::DFT_one()
{
m_dft_vector = NULL;
m_has_dft_vector = false;
m_dft_vector_size = 0;
}
DFT_one::~DFT_one(void)
{
if (m_has_dft_vector && (NULL != m_dft_vector) && (m_dft_vector_size>0))
delete[] m_dft_vector;
}
bool DFT_one::has_dft_vector()
{
return m_has_dft_vector;
}
void DFT_one::clear_dft_vector()
{
if (m_has_dft_vector && (NULL != m_dft_vector) && (m_dft_vector_size>0)) {
delete[] m_dft_vector;
m_has_dft_vector = false;
m_dft_vector_size = 0;
}
}
void DFT_one::clear_idft_vector()
{
if (m_has_idft_vector && (NULL != m_idft_vector) && (m_idft_vector_size>0)) {
delete[] m_idft_vector;
m_has_idft_vector = false;
m_idft_vector_size = 0;
}
}
void DFT_one::print()
{
if ((!m_has_dft_vector) || (NULL == m_dft_vector) || (m_dft_vector_size <= 0))
return;
for (int i = 0; i qDebug()< }
for (int i = 0; i qDebug()< }
}
bool DFT_one::fft(Complex vec[], int len)
{
this->clear_dft_vector();
this->d1_vector.clear();
int old_len=len,r=0;
if((len<=0) || (NULL==vec))
return false;
clear_dft_vector();
d1_vector.clear();
if(len&(len-1)==0){
int k=1;
while(k!=len){
r++;
k=k*2;
}
}else{
int k=1;
while(k r++;
k=k*2;
}
len=k;
}
m_dft_vector = new Complex[len];
Complex *cp=new Complex[len];
for (int u = 0; u m_dft_vector[u].SetValue(0, 0);
if(u cp[u]=vec[u];
}else{
cp[u]=Complex(0,0);
}
}
int i,j,k; // 循环变量
int bfsize,p;
double angle; // 角度
Complex *W,*X1,*X2,*X;
len = 1 < // 分配运算所需存储器
W = new Complex[len / 2];
X1 = new Complex[len];
X2 = new Complex[len];
// 计算加权系数
for(i = 0; i {
angle = -i * PI * 2 / len;
W[i] = Complex (cos(angle), sin(angle));
}
// 将时域点写入X1
for(int i=0;i X1[i]=cp[i];
}
// 采用蝶形算法进行快速付立叶变换
for(k = 0; k {
for(j = 0; j <1 < {
bfsize = 1 <<(r-k);//做蝶形运算两点间距离
for(i = 0; i {
p = j * bfsize;
X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2])
* W[i * (1< }
}
X = X1;
X1 = X2;
X2 = X;
}
// 重新排序
for(j = 0; j {
p = 0;
for(i = 0; i {
if (j&(1< {
p+=1<<(r-i-1);
}
}
m_dft_vector[j]=X1[p];
d1_vector.push_back(X1[p]);
}
delete W;
delete X1;
delete X2;
m_has_dft_vector = true;
m_dft_vector_size = len;
return true;
}
bool DFT_one::ifft(Complex vec[], int len)
{
int old_len=len,r=0;
if((len<=0) || (NULL==vec))
return false;
clear_idft_vector();
id1_vector.clear();
if(len&(len-1)==0){
int k=1;
while(k!=len){
r++;
k=k*2;
}
}else{
int k=1;
while(k r++;
k=k*2;
}
len=k;
}
m_idft_vector = new Complex[len];
Complex *cp=new Complex[len];
for (int u = 0; u m_idft_vector[u].SetValue(0, 0);
if(u cp[u]=vec[u];
}else{
cp[u]=Complex(0,0);
}
}
int i,j,k; // 循环变量
int bfsize,p;
double angle; // 角度
Complex *W,*X1,*X2,*X;
len = 1 < // 分配运算所需存储器
W = new Complex[len / 2];
X1 = new Complex[len];
X2 = new Complex[len];
// 计算加权系数
for(i = 0; i {
angle = -i * PI * 2 / len;
W[i] = Complex (cos(angle), -sin(angle));
}
// 将时域点写入X1
for(int i=0;i X1[i]=cp[i];
}
// 采用蝶形算法进行快速付立叶变换
for(k = 0; k {
for(j = 0; j <1 < {
bfsize = 1 <<(r-k);//做蝶形运算两点间距离
for(i = 0; i {
p = j * bfsize;
X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2])
* W[i * (1< }
}
X = X1;
X1 = X2;
X2 = X;
}
// 重新排序
for(j = 0; j {
p = 0;
for(i = 0; i {
if (j&(1< {
p+=1<<(r-i-1);
}
}
m_idft_vector[j]=X1[p]/len;
id1_vector.push_back(X1[p]/len);
}
delete W;
delete X1;
delete X2;
m_has_idft_vector = true;
m_idft_vector_size = len;
return true;
}
bool DFT_one::fft(double vec[], int len)
{
Complex *v=new Complex[len];
for(int i=0;i v[i]=Complex(vec[i],0);
}
bool fd=fft(v,len);
return fd;
}


dft_two.h

#ifndef DFT_TWO_H
#define DFT_TWO_H
#include "complex.h"
#include "dft_one.h"
#include
#include
#include
#include
class DFT_two
{
public:
DFT_two(void);
~DFT_two(void);
public:
bool fft2(double matrix[], int width, int height);
bool ifft2(Complex matrix[], int width, int height);
void re_normalize_spectrum(double fp[], double width, double height);
bool has_dft2_matrix(); // 是否已存有变换结果
void clear_dft2_matrix(); // 清除已有的变换结果
void clear_idft2_matrix();
void print_matrix(); // 打印变换结果
public:
Complex *m_dft2_matrix;
double max;
QVector d2_matrix;
QVector d2_data;
Complex *m_idft2_matrix;
QVector id2_data;
double *re_m_spectrum_data;
protected:
bool m_has_dft_matrix;
bool m_is_normalized;
bool m_is_spectrum_shifted;
int m_dft_matrix_height;
int m_dft_matrix_width;
bool m_has_idft_matrix;
int m_idft_matrix_height;
int m_idft_matrix_width;
};
#endif // DFT_TWO_H


dft_two.cpp

#include "dft_two.h"
DFT_two::DFT_two(){
m_dft2_matrix = NULL;
m_has_dft_matrix = false;
m_is_normalized = false;
m_is_spectrum_shifted = false;
m_dft_matrix_height = 0;
m_dft_matrix_width = 0;
}
DFT_two::~DFT_two(void)
{
if (m_has_dft_matrix && (NULL != m_dft2_matrix) && ((m_dft_matrix_height*m_dft_matrix_width)>0))
delete[] m_dft2_matrix;
}
bool DFT_two::has_dft2_matrix()
{
return m_has_dft_matrix;
}
void DFT_two::clear_dft2_matrix()
{
if (m_has_dft_matrix && (NULL != m_dft2_matrix) && ((m_dft_matrix_height*m_dft_matrix_width)>0)) {
delete[] m_dft2_matrix;
m_has_dft_matrix = false;
m_dft_matrix_height = 0;
m_dft_matrix_width = 0;
}
}
void DFT_two::clear_idft2_matrix()
{
if (m_has_idft_matrix && (NULL != m_idft2_matrix) && ((m_idft_matrix_height*m_idft_matrix_width)>0)) {
delete[] m_idft2_matrix;
m_has_idft_matrix = false;
m_idft_matrix_height = 0;
m_idft_matrix_width = 0;
}
}
void DFT_two::print_matrix()
{
if ((!m_has_dft_matrix) || (NULL == m_dft2_matrix) || (m_dft_matrix_height <= 0) || (m_dft_matrix_width <= 0))
return;
for (int u = 0; u for (int v = 0; v qDebug()< }
}
}
bool DFT_two::fft2( double matrix[], int width, int height){
if (((width*height) <= 0) || (NULL == matrix))
return false;
clear_dft2_matrix();
d2_matrix.clear();
m_dft2_matrix = new Complex[width*height];
DFT_one d1;
for(int i=0;i d1.clear_dft_vector();
d1.fft(matrix+i*width,width);
for(int j=0;j m_dft2_matrix[j+i*width]=d1.d1_vector[j];
d2_matrix.push_back(d1.d1_vector[j]);
}
}
for(int i=0;i Complex *yl=new Complex[height];
for(int j=0;j yl[j]=m_dft2_matrix[i+j*width];
}
d1.clear_dft_vector();
d1.fft(yl,height);
for(int j=0;j m_dft2_matrix[i+j*width]=d1.d1_vector[j];
d2_matrix[i+j*width]=d1.d1_vector[j];
}
}
m_has_dft_matrix = true;
m_dft_matrix_height = height;
m_dft_matrix_width = width;
return true;
}
bool DFT_two::ifft2( Complex matrix[], int width, int height){
if (((width*height) <= 0) || (NULL == matrix))
return false;
clear_idft2_matrix();
m_idft2_matrix = new Complex[width*height];
DFT_one d1;
for(int i=0;i for(int j=0;j m_idft2_matrix[j+i*width]=matrix[i*width+j];
}
}
for(int i=0;i Complex *yl=new Complex[height];
for(int j=0;j yl[j]=m_idft2_matrix[i+j*width];
}
d1.clear_idft_vector();
d1.ifft(yl,height);
for(int j=0;j m_idft2_matrix[i+j*width]=d1.id1_vector[j];
}
}
for(int i=0;i d1.clear_idft_vector();
d1.ifft(m_idft2_matrix+i*width,width);
for(int j=0;j m_idft2_matrix[j+i*width]=d1.id1_vector[j];
}
}
m_has_idft_matrix = true;
m_idft_matrix_height = height;
m_idft_matrix_width = width;
return true;
}


拿Qt做的一个图像增强

mainWindow.h

#ifndef MAINWINDOW_H
#define MAINWINDOW_H
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include "complex.h"
#include "dft_one.h"
#include "dft_two.h"
namespace Ui {
class MainWindow;
}
class MainWindow : public QMainWindow
{
Q_OBJECT
public:
explicit MainWindow(QWidget *parent = 0);
~MainWindow();
private slots:
void on_choiceButton_clicked();
void on_paintButton_clicked();
void on_paintButton_2_clicked();
private:
Ui::MainWindow *ui;
QString path;
QVectorimg;
int width;
int height;
QImage *res;
};
#endif // MAINWINDOW_H


mainWindow.cpp

#include "mainwindow.h"
#include "ui_mainwindow.h"
MainWindow::MainWindow(QWidget *parent) :
QMainWindow(parent),
ui(new Ui::MainWindow)
{
ui->setupUi(this);
}
MainWindow::~MainWindow()
{
delete ui;
}
void MainWindow::on_choiceButton_clicked()
{
path = QFileDialog::getOpenFileName(this, tr("open image file"),
"./", tr("Image files(*.bmp *.jpg *.pbm *.pgm *.png *.ppm *.xbm *.xpm);;All files (*.*)"));
QPixmap p(path);
QImage image = p.toImage();
if(path.isEmpty()==false)
{
img.clear();
// ui->imgLabel->setPixmap(p.scaled(ui->imgLabel->width(),
// ui->imgLabel->height(), Qt::KeepAspectRatio, Qt::SmoothTransformation));
ui->imgLabel->setPixmap(p);
警告","打开图片失败!");
return;
}
}
double H(double D2,double rH,double rL,int D0){
double c=1;
return (rH-rL)*(1-exp(-c*D2/pow(D0,2)))+rL;
}
void MainWindow::on_paintButton_clicked()
{
double rH=5,rL=0.95;
int D0=200;
int xk=2,yk=2;
while(1){
if(pow(2,xk-1)=width) break;
xk++;
}
while(1){
if(pow(2,yk-1)=height) break;
yk++;
}
int new ';
d2.ifft2(d2.m_dft2_matrix,newWidth,newHeight);
for(int i=0;i if(d2.m_idft2_matrix[i].real>255) d2.m_idft2_matrix[i].real=255;
if(d2.m_idft2_matrix[i].real<0) d2.m_idft2_matrix[i].real=0;
if(matrix[i]>255) matrix[i]=255;
if(matrix[i]<0) matrix[i]=0;
}
QImage *hz=new QImage(width, height, QImage::Format_ARGB32);
for(int i=0;i for(int j=0;j hz->setPixel(j,i,qRgb(d2.m_idft2_matrix[j+i*newWidth].real,
d2.m_idft2_matrix[j+i*newWidth].real,
d2.m_idft2_matrix[j+i*newWidth].real));
}
}
// ui->imgLabel_2->setPixmap(QPixmap::fromImage(*hz).scaled(ui->imgLabel_2->width(),
// ui->imgLabel_2->height(), Qt::KeepAspectRatio, Qt::SmoothTransformation));
ui->imgLabel_2->setPixmap(QPixmap::fromImage(*hz));
//用于保存
res=hz;
}
void MainWindow::on_paintButton_2_clicked()
{
QString filename1 = QFileDialog::getSaveFileName(this,tr("Save Image"),"图片.bmp",tr("Images (*.png *.bmp *.jpg)")); //选择路径
res->save(filename1);
}


原图

处理后