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

处理后

C++ 实现二维fft和ifft

上一篇:旋转矩阵----列表顺时针旋转90°


下一篇:在行列都排好序的矩阵中找指定的数