6.3.2巴特沃斯（butterworth）低通滤波器
在本程序中，设计了巴特沃斯滤频域波器函数，可以实现频域的低通、高通、带通和带阻滤波。函数如下：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};
在本程序中，共有六个自定义函数，分别是：
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=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时的滤波结果：
版权声明：本文为phoenixdsg原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。