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 <QDebug>
#include <QString>
#include <qmath.h>
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 <QVector>
#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;
QVector<Complex>d1_vector;
QVector<Complex>id1_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<m_dft_vector_size; i++) {
qDebug()<<m_dft_vector[i]<<" ";
}
for (int i = 0; i<m_idft_vector_size; i++) {
qDebug()<<m_idft_vector[i]<<" ";
}
}
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<len){
r++;
k=k*2;
}
len=k;
}
m_dft_vector = new Complex[len];
Complex *cp=new Complex[len];
for (int u = 0; u < len; u++) {
m_dft_vector[u].SetValue(0, 0);
if(u<old_len){
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 << r; // 计算付立叶变换点数
// 分配运算所需存储器
W = new Complex[len / 2];
X1 = new Complex[len];
X2 = new Complex[len];
// 计算加权系数
for(i = 0; i < len / 2; i++)
{
angle = -i * PI * 2 / len;
W[i] = Complex (cos(angle), sin(angle));
}
// 将时域点写入X1
for(int i=0;i<len;i++){
X1[i]=cp[i];
}
// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < r; k++)//k为蝶形运算的级数
{
for(j = 0; j < 1 << k; j++)
{
bfsize = 1 << (r-k);//做蝶形运算两点间距离
for(i = 0; i < bfsize / 2; 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<<k)];
}
}
X = X1;
X1 = X2;
X2 = X;
}
// 重新排序
for(j = 0; j < len; j++)
{
p = 0;
for(i = 0; i < r; i++)
{
if (j&(1<<i))
{
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<len){
r++;
k=k*2;
}
len=k;
}
m_idft_vector = new Complex[len];
Complex *cp=new Complex[len];
for (int u = 0; u < len; u++) {
m_idft_vector[u].SetValue(0, 0);
if(u<old_len){
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 << r; // 计算付立叶变换点数
// 分配运算所需存储器
W = new Complex[len / 2];
X1 = new Complex[len];
X2 = new Complex[len];
// 计算加权系数
for(i = 0; i < len / 2; i++)
{
angle = -i * PI * 2 / len;
W[i] = Complex (cos(angle), -sin(angle));
}
// 将时域点写入X1
for(int i=0;i<len;i++){
X1[i]=cp[i];
}
// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < r; k++)//k为蝶形运算的级数
{
for(j = 0; j < 1 << k; j++)
{
bfsize = 1 << (r-k);//做蝶形运算两点间距离
for(i = 0; i < bfsize / 2; 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<<k)];
}
}
X = X1;
X1 = X2;
X2 = X;
}
// 重新排序
for(j = 0; j < len; j++)
{
p = 0;
for(i = 0; i < r; i++)
{
if (j&(1<<i))
{
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<len;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 <QDebug>
#include <QVector>
#include <cmath>
#include <cstring>
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<Complex> d2_matrix;
QVector<double> d2_data;
Complex *m_idft2_matrix;
QVector<double> 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<m_dft_matrix_height; u++) {
for (int v = 0; v<m_dft_matrix_width; v++) {
qDebug()<<m_dft2_matrix[v + u*m_dft_matrix_width];
}
}
}
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<height;i++){
d1.clear_dft_vector();
d1.fft(matrix+i*width,width);
for(int j=0;j<width;j++){
m_dft2_matrix[j+i*width]=d1.d1_vector[j];
d2_matrix.push_back(d1.d1_vector[j]);
}
}
for(int i=0;i<width;i++){
Complex *yl=new Complex[height];
for(int j=0;j<height;j++){
yl[j]=m_dft2_matrix[i+j*width];
}
d1.clear_dft_vector();
d1.fft(yl,height);
for(int j=0;j<height;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<height;i++){
for(int j=0;j<width;j++){
m_idft2_matrix[j+i*width]=matrix[i*width+j];
}
}
for(int i=0;i<width;i++){
Complex *yl=new Complex[height];
for(int j=0;j<height;j++){
yl[j]=m_idft2_matrix[i+j*width];
}
d1.clear_idft_vector();
d1.ifft(yl,height);
for(int j=0;j<height;j++){
m_idft2_matrix[i+j*width]=d1.id1_vector[j];
}
}
for(int i=0;i<height;i++){
d1.clear_idft_vector();
d1.ifft(m_idft2_matrix+i*width,width);
for(int j=0;j<width;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 <QMainWindow>
#include <QWidget>
#include <QFile>
#include <QDebug>
#include <QPushButton>
#include <QDateTime>
#include <QString>
#include <QMessageBox>
#include <QFileDialog>
#include <QFileInfo>
#include <QLabel>
#include <QVector>
#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;
QVector<double>img;
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);
width=image.width();
height=image.height();
for(int i=0;i<image.height();i++){
for(int j=0;j<image.width();j++){
QColor Colos= image.pixelColor(j,i);
int degree=qGray(Colos.red(),Colos.green(),Colos.blue());//单通道灰度
img.push_back((double)degree);
}
}
}else{
QMessageBox mesg;
mesg.warning(this,"警告","打开图片失败!");
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&&pow(2,xk)>=width) break;
xk++;
}
while(1){
if(pow(2,yk-1)<height&&pow(2,yk)>=height) break;
yk++;
}
int newWidth=pow(2,xk);
int newHeight=pow(2,yk);
//对图片进行补零操作
double *matrix=new double[newWidth*newHeight];
int imgIndex=0;
for(int i=0;i<newWidth*newHeight;i++){
if(i%newWidth<width&&(i/newWidth)<height){
matrix[i]=img.at(imgIndex);
imgIndex++;
}else{
matrix[i]=0;
}
}
DFT_two d2;
d2.fft2(matrix,newWidth,newHeight);
for(int i=0;i<newHeight;i++){
for(int j=0;j<newWidth;j++){
double D2=0;
if(i<newHeight/2&&j<newWidth/2){
D2=pow(i,2)+pow(j,2);
}else if(i>=newHeight/2&&j<newWidth/2){
D2=pow(newHeight-i,2)+pow(j,2);
}else if(i<newHeight/2&&j>=newWidth/2){
D2=pow(i,2)+pow(newWidth-j,2);
}else if(i>=newHeight/2&&j>=newWidth/2){
D2=pow(newHeight-i,2)+pow(newWidth-j,2);
}
d2.m_dft2_matrix[j+i*newWidth]=d2.m_dft2_matrix[j+i*newWidth]*Complex(H(D2,rH,rL,D0),0);
}
}
qDebug()<<' ';
d2.ifft2(d2.m_dft2_matrix,newWidth,newHeight);
for(int i=0;i<newWidth*newHeight;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<height;i++){
for(int j=0;j<width;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);
}
原图
处理后
标签：ifft,matrix,dft,fft,C++,vector,Complex,DFT,include 来源： https://www.cnblogs.com/God-Destiny/p/15744548.html