C/C++教程
C++ 实现二维fft和ifft
本文主要是介绍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);
}
原图
处理后
这篇关于C++ 实现二维fft和ifft的文章就介绍到这儿，希望我们推荐的文章对大家有所帮助，也希望大家多多支持为之网！
您可能喜欢
-
怎么实现ansible playbook 备份代码中命名包含时间戳功能？-icode9专业技术文章分享11-22
-
ansible 的archive 参数是什么意思？-icode9专业技术文章分享11-22
-
ansible 中怎么只用archive 排除某个目录？-icode9专业技术文章分享11-22
-
exclude_path参数是什么作用？-icode9专业技术文章分享11-22
-
微信开放平台第三方平台什么时候调用数据预拉取和数据周期性更新接口？-icode9专业技术文章分享11-22
-
uniapp 实现聊天消息会话的列表功能怎么实现？-icode9专业技术文章分享11-22
-
在Mac系统上将图片中的文字提取出来有哪些方法？-icode9专业技术文章分享11-22
-
excel 表格中怎么固定一行显示不滚动？-icode9专业技术文章分享11-22
-
怎么将 -rwxr-xr-x 修改为 drwxr-xr-x？-icode9专业技术文章分享11-22
-
在Excel中怎么将小数向上取整到最接近的整数？-icode9专业技术文章分享11-22
-
Excel中常见的取整函数有哪些？-icode9专业技术文章分享11-22
-
错误信息：does not match the actual quantity supplied 提示是什么意思？-icode9专业技术文章分享11-22
-
WebSocket学习：初学者的简单指南11-21
-
获取apk的md5值有哪些方法？-icode9专业技术文章分享11-20
-
xml报文没有传 IdentCode ，为什么正常解析没报错呢？-icode9专业技术文章分享11-20
栏目导航
