首页 > 代码库 > 6.3.2巴特沃斯(butterworth)低通滤波器

6.3.2巴特沃斯(butterworth)低通滤波器

 

/**************程序说明****************
在本程序中,共有六个自定义函数,分别是:
1. myMagnitude(Mat & complexImg,Mat & magnitudeImage),在该函数中封装了Opencv中的
   magnitude函数,实现对于复数图像的幅值计算。该函数共有两个参数:
   complexImg--输入的复数阵列,或复数图像
   magnitudeImage--输出的幅值阵列,或幅值图像

2. dftshift(Mat& ds),该函数实现对图像四个象限的对角互换,相当于MatLab中 fftshift(),将频谱的原点(0,0)移到图像中心。
3. srcCentralized(Mat& src)用于傅里叶变换前的预处理,以便得到傅里叶频谱的原点(0,0)位于图像的中心。
   该函数与dftshift()目的一致,实现方法不同,一个是变换前预处理,一个是变换后处理。
4. imshowComplexMat(Mat&dftDst,String winName,bool inverseSpectrum),该函数用于显示复数图像或双通道矩阵,共有三个参数:
   dftDst--待显示的复数矩阵
   winName--显示复数矩阵的窗口名字
   inverseSpectrum-输入的dftDst是正向傅里叶变换的结果,还是逆傅里叶变换的结果
5. createFilterButterworth(Mat&filter,int n,int R,int W,FilterForm filterform),用于制作Butterworth频域滤波器,该函数利用了ptr()
   指针遍历图像的方法。该函数可以实现低通、高通、带通、带阻滤波器。目前该函数共有五个参数:
   filter--输入的矩阵,要求数据类型为CV_64FC2;
   n--巴特沃斯阶数
   R--截止频率半径,如果小于0,则返回一个全口径滤波器,否则返回一个口径受限的滤波器
   W--带宽
   filterform--滤波器形式,它是个枚举类型数据,enum FilterForm{LOW_PASS_FILTER,HIGH_PASS_FILTER,BAND_PASS_FILTER,BAND_STOP_FILTER};

6.void myDft(Mat&src,Mat&dst,bool isProCentralized,bool doubleSizeOrNot),该函数更有四个参数
  src--是输入的原图像
  dst--是傅里叶变换的输出图像
  isProCentralized--表示是否调用SRCCentralized函数,对src进行中心化预处理
  doubleSizeOrNot--表示是否需要将原图像尺寸扩展为两倍,以便解决卷积缠绕问题

 

 

#include <iostream>
#include<opencv2/opencv.hpp>
using namespace cv;
using namespace std;
#define CV_MAT_ELEM2(src,dtype,y,x) \
    (dtype*)(src.data+y*src.step[0]+x*src.step[1])

/**************程序说明****************

在主程序内容,参见注释。
**********************************************/
/*1.求取复数矩阵的幅值*/

void myMagnitude(Mat & complexImg,Mat & magnitudeImage)
{
    Mat planes[2];
    split(complexImg,planes);
    magnitude(planes[0],planes[1],magnitudeImage);
}

/*傅里叶变换后的频谱图后处理,将傅里叶普的原点(0,0)平移到图像的中心*/
void dftshift(Mat& ds)
{
    int cx=ds.cols/2;//图像的中心点x 坐标
    int cy=ds.rows/2;//图像的中心点y 坐标
    Mat q0=ds(Rect(0,0,cx,cy));//左上
    Mat q1=ds(Rect(cx,0,cx,cy));//右上
    Mat q2=ds(Rect(0,cy,cx,cy));//左下
    Mat q3=ds(Rect(cx,cy,cx,cy));//右下
    Mat tmp;
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);
    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);
}
/*傅里叶变换前的预处理,以便频谱图的原点(0,0)移动到图像的中心*/
void srcCentralized(Mat& src)
{
    switch (src.depth())
    {
    case CV_32F:
        for(int i=0;i<src.rows;i++)
        {
            for(int j=0;j<src.cols;j++)
            {
                float* mv=CV_MAT_ELEM2(src,float,i,j);
                if((i+j)%2!=0)
                {
                    for(int c=0;c<src.channels();c++)//对所有通道同样操作
                        mv[c]=-mv[c];//如果i+j为奇数,该像素值取负值
                }
            }
        }
        break;
    case CV_64F:
        for(int i=0;i<src.rows;i++)
        {
            for(int j=0;j<src.cols;j++)
            {
                double* mv=CV_MAT_ELEM2(src,double,i,j);
                if((i+j)%2!=0)
                {
                    for(int c=0;c<src.channels();c++)//遍历各个通道
                        mv[c]=-mv[c];//如果i+j为奇数,该像素值取负值
                }
            }
        }
        break;

    default:
        break;
    }

}
/*在窗口中显示复数图像,如果是正向傅里叶矩阵,需要取log才能显示更多频谱信息
 *如果是逆傅里叶变换,通过normalize归一化后,显示频谱图*/
void imshowComplexMat(Mat&dftDst,String winName,bool inverseSpectrum)
{
    Mat magI;
    myMagnitude(dftDst,magI);
    if(!inverseSpectrum)//如果是正向傅里叶变换谱
    {
        magI+=Scalar::all(1);
        log(magI,magI);
    }
    normalize(magI,magI,0,1,NORM_MINMAX);
    imshow(winName,magI);
}

///
enum FilterForm{LOW_PASS_FILTER,HIGH_PASS_FILTER,BAND_PASS_FILTER,BAND_STOP_FILTER};
void createFilterButterworth(Mat&filter,int n,int R,int W,FilterForm filterform)
{
    double Rs=R*R;//R1_square
    int cx=filter.cols/2;
    int cy=filter.rows/2;
    switch(filterform)
    {
    case LOW_PASS_FILTER:
        for(int i=0;i<filter.rows;i++)
        {
            Vec2d* pf=filter.ptr<Vec2d>(i);
            for(int j=0;j<filter.cols;j++)
            {
                double rs=(j-cx)*(j-cx)+(i-cy)*(i-cy);//rs表示r的平方
                pf[j][0]=1./(1.+pow(rs/Rs,n));//Rs是R的平方,
                pf[j][1]=pf[j][0];
            }
        }
        break;
    case HIGH_PASS_FILTER:
        for(int i=0;i<filter.rows;i++)
        {
            Vec2d* pf=filter.ptr<Vec2d>(i);
            for(int j=0;j<filter.cols;j++)
            {
                double rs=(j-cx)*(j-cx)+(i-cy)*(i-cy);
                double Lp=1./(1.+pow(rs/Rs,n));//巴特沃斯公式
                pf[j][0]=1.0-Lp;
                pf[j][1]=pf[j][0];
            }
        }
        break;
    case BAND_STOP_FILTER:
        for(int i=0;i<filter.rows;i++)
        {
            Vec2d* pf=filter.ptr<Vec2d>(i);
            for(int j=0;j<filter.cols;j++)
            {
                double rs=(j-cx)*(j-cx)+(i-cy)*(i-cy);
                double r=std::sqrt(rs);//r相当于书上的D
                pf[j][0]=1./(1.+pow(r*W/(rs-Rs),2*n));
                pf[j][1]=pf[j][0];
            }
        }
        break;
    case BAND_PASS_FILTER:
        for(int i=0;i<filter.rows;i++)
        {
            Vec2d* pf=filter.ptr<Vec2d>(i);
            for(int j=0;j<filter.cols;j++)
            {
                double rs=(j-cx)*(j-cx)+(i-cy)*(i-cy);
                double r=std::sqrt(rs);//r相当于书上的D
                pf[j][0]=1-1./(1.+pow(r*W/(rs-Rs),2*n));
                pf[j][1]=pf[j][0];
            }
        }
        break;
    }
///***************【显示滤波器,如果不需要显示,将代码注销】*************/
    Mat displayFilter;
    extractChannel(filter,displayFilter,1);
    imshow("btw filter image",displayFilter);

}


void myDft(Mat&src,Mat&dst,bool isProCentralized=false,bool doubleSizeOrNot=false)
{
    CV_Assert(src.channels()==1);//验证src是否是单通道
    int ny=src.rows,nx=src.cols;
    Mat srcPadded;
    if(doubleSizeOrNot)//如果doubleOrNot为真,对src尺度扩展一倍
    {
        cv::copyMakeBorder(src,srcPadded,0,ny,0,nx,BORDER_CONSTANT);
    }
    else//否则,将src填补为最优傅里叶变换尺寸
    {
        int padx=getOptimalDFTSize(nx);//获取最优傅里叶变换列数
        int pady=getOptimalDFTSize(ny);//获取最优傅里叶变换行数数
        cv::copyMakeBorder(src,srcPadded,0,pady-ny,0,padx-nx,BORDER_CONSTANT);
    }
    if(srcPadded.type()!=CV_64FC1)
        srcPadded.convertTo(srcPadded,CV_64FC1);//转成浮点数据类型
    //如果isProCentralized为真,则在傅里叶变换前,对src做中心化预处理,
    //这样在傅里叶变换后的频谱图的(0,0)点就会位于频谱图的中心
    if(isProCentralized)
        srcCentralized(srcPadded);

    Mat planes[2]={srcPadded,Mat::zeros(srcPadded.rows,srcPadded.cols,srcPadded.type())};
    Mat srcComplex;
    merge(planes,2,srcComplex);
    dft(srcComplex,dst,DFT_COMPLEX_OUTPUT,0);//离散傅立叶变换
    cout<<"srcPadded.rows="<<srcPadded.rows<<"  "<<"src.cols="<<srcPadded.cols<<endl;
}

int main()
{
    ///读入灰度图像
    Mat src=http://www.mamicode.com/imread("D:\\Qt\\MyImage\\Fig0333(a).tif",0);

    ///***************【调用自定义的傅里叶变换myDft()】*************/
    Mat dftDst;//预声明dft的输出结果矩阵
    myDft(src,dftDst,0,0);//调用自定义的dft函数,对src执行傅里叶变换,dftDst是傅里叶变换结果
    dftshift(dftDst);//将傅里叶变换的结果,四象限对角互换
    imshowComplexMat(dftDst,"dftSpectrum display",false);//显示傅里叶频谱图

    ///***************【创建滤波器】*************/
    Mat filter(dftDst.rows,dftDst.cols,CV_64FC2,Scalar(0,0));
//巴特沃斯滤波器的阶数n=2,半径分别取R=10,30,60,160,460 createFilterButterworth(filter,
2,100,60,LOW_PASS_FILTER); ///***************【频域滤波】*************/ Mat temp; temp=dftDst.mul(filter); ///***************【调用opencv中的dft,执行逆傅里叶变换】*************/ Mat filteredResult; dft(temp,filteredResult,DFT_INVERSE); imshowComplexMat(filteredResult,"Frequency domain butterworth filtered image",1);
imshow("srcimage",src);
    waitKey();
    return 0;
}

 技术分享      技术分享

技术分享    技术分享

 上面四幅图,分别是巴特沃斯滤波器,n=2,截止半径=10;第二副图像是src的频谱图;第三幅图像是src原图像;第四幅是滤波后的结果。

 下面是截止半径R分别为30,60,160和460时的滤波结果:

技术分享   技术分享

技术分享    技术分享

 

 

 

 

 

 

 

 

 

 

 

 

 


6.3.2巴特沃斯(butterworth)低通滤波器