理想低通滤波器

     理想低通滤波器：该系统如下：

 

H(u,v)=1  D(u,v)<=D₀

H(u,v)=0  D(u,v)>D₀

其中D(u,v)=[(u-M/2)² +(v-n/2)² ]^1/2

     对一幅图像的处理，在岗萨雷斯的数字图像处理中有说明，在这里给帖出来：

一，用(-1)^(x+y) 乘以输入图像来进行中心变换；

二，由(1)计算图像的DFT，即F(u,v);

三，用滤波器函数H(u,v)乘以F(u,v)；

四，计算(3)中结果的反DFT;

五，得到(4)中结果的实部；

六，用(-1)^(x+y) 乘以(5)中的结果；

这样就把理想低通滤波器实现了；

 

关于理想低通滤波器的时域表示，下图列出了不同截止频率时系统的空间图像以及使用相应的截止频率进行滤波时对应的结果：

         

截止频率为5          对应滤波结果

           

截止频率为15        对应滤波结果 

          

截止频率为30        对应滤波结果  

           

截止频率为80        对应滤波结果 

            

截止频率为230       对应滤波结果

 

      可以明显的看出，使用理想低通滤波，结果图像会产生较大的振铃，而且低理想滤波器通过物理器件是实现不了的；

      关于频率滤波，可以使用两种方法来实现；上面已经给出了一种；另一种方法是把滤波器的系统进行反傅里叶变换，然后将结果与图像进行卷积，关于这一块，我测试了一下，目前来说还存在一些问题，过些时日我把这一块再写出来!

      下面把代码放帖出来：

 

.h文件 

 // IdealLowPass.h: interface for the CIdealLowPass class. // // #if !defined(AFX_IDEALLOWPASS_H__3E7EC990_C6A1_46EF_8051_26586B835D5F__INCLUDED_) #define AFX_IDEALLOWPASS_H__3E7EC990_C6A1_46EF_8051_26586B835D5F__INCLUDED_ #if _MSC_VER > 1000 #pragma once #endif // _MSC_VER > 1000 /********************************************** *类名:理想低通滤波器 *目的:对一幅图像进行低通滤波 *作者:无忧狂澜 *日期:2010-02-06 **********************************************/ //低通滤波器要用到卷积和FFT #include "Convolution.h" #include "FFT.h" #define X 256 #define Y 256 #define U 256 #define V 256 #define PI 3.1415926 #define TWOPI 6.2831852 class CIdealLowPass { public: void ComputeSpaceSystem(); CIdealLowPass(); virtual ~CIdealLowPass(); private: double m_fHXY_Real[X*Y]; //空间系统复数 double m_fHXY_Imag[X*Y]; double m_fHUV_Real[U*V]; //频域系统复数 double m_fHUV_Imag[U*V]; double m_fHUVTotal; //u,v里的总和,为了计算方便起见才做的 int m_nMSampling; //采样点的个数 int m_nNSampling; //采样点的个数 double m_fCutFreq; //截止频率D0 CFFT m_fft; //傅里叶变换及逆变换 public: void SetMSampling(int nMSampling) {m_nMSampling=nMSampling;} void SetNSampling(int nNSampling) {m_nNSampling=nNSampling;} void SetCutFreq(double fCutFreq) {m_fCutFreq=fCutFreq;} double GetHXYReal(int nX=0,int nY=0) {return m_fHXY_Real[nX*Y+nY];} double GetHXYImag(int nX=0,int nY=0) {return m_fHXY_Imag[nX*Y+nY];} /*****************不使用卷积滤波方法begin*******************/ void LowPassFilter(LPSTR lpDIBBits, LONG lLineBytes, LONG lWidth, LONG lHeight); //低通滤波 /*****************不使用卷积滤波方法end*******************/ }; #endif // !defined(AFX_IDEALLOWPASS_H__3E7EC990_C6A1_46EF_8051_26586B835D5F__INCLUDED_)

 

.cpp文件

// IdealLowPass.cpp: implementation of the CIdealLowPass class. // // #include "stdafx.h" #include "mydip.h" #include "IdealLowPass.h" #include <math.h> #ifdef _DEBUG #undef THIS_FILE static char THIS_FILE[]=__FILE__; #define new DEBUG_NEW #endif // // Construction/Destruction // CIdealLowPass::CIdealLowPass() { //m_nMSampling=nMSampling; //采样点的个数 //m_nNSampling=nNSampling; //采样点的个数 //m_fCutFreq=fCutFreq; //截止频率D0 } CIdealLowPass::~CIdealLowPass() { } //计算空间系统 void CIdealLowPass::ComputeSpaceSystem() { //先计算频域系统m_fHUV double fTemp=0; m_fHUVTotal=0; for(int u=0;u<m_nMSampling;u++) for(int v=0;v<m_nNSampling;v++) { fTemp=sqrt((u-m_nMSampling/2)*(u-m_nMSampling/2)+(v-m_nNSampling/2)*(v-m_nNSampling/2)); if(fTemp<=m_fCutFreq) m_fHUV_Real[v+m_nNSampling*u]=1; else m_fHUV_Real[v+m_nNSampling*u]=0; //m_fHUVTotal+=m_fHUV[u][v]; } //计算空间域系统m_fHXY //反向FFT m_fft.ComputeConverseFFT(m_nMSampling,m_nNSampling,m_fHUV_Real,NULL,m_fHXY_Real,m_fHXY_Imag,false); //得到空间域数据 /* for(int x=0;x<X;x++) for(int y=0;y<Y;y++) { m_fHXY[x][y]=0; for(int u=0;u<m_nMSampling;u++) for(int v=0;v<m_nNSampling;v++) { m_fHXY[x][y]+=m_fHUV[u][v]*tan(2*PI*(u*x/m_nMSampling+v*y/m_nNSampling)); } } */ } void CIdealLowPass::LowPassFilter(LPSTR lpDIBBits, LONG lLineBytes, LONG lWidth, LONG lHeight) { double *fFReal,*fFImag,*fTReal,*fTImag; double fTemp=0; int i,j; unsigned char* lpSrc; fFReal=new double[lWidth*lHeight]; fFImag=new double[lWidth*lHeight]; fTReal=new double[lWidth*lHeight]; fTImag=new double[lWidth*lHeight]; //1,中心化图像FFT m_fft.Fourier(lpDIBBits,lWidth,lHeight,fFReal,fFImag,true); //2,F*H //计算H(u,v) for(int u=0;u<lWidth;u++) for(int v=0;v<lHeight;v++) { fTemp=sqrt((u-m_nMSampling/2)*(u-m_nMSampling/2)+(v-m_nNSampling/2)*(v-m_nNSampling/2)); if(fTemp<=m_fCutFreq) m_fHUV_Real[v+m_nNSampling*u]=1; else m_fHUV_Real[v+m_nNSampling*u]=0; } //计算F*H for(i=0;i<lHeight;i++) for(j=0;j<lWidth;j++) { fFReal[j*lHeight+i]=fFReal[j*lHeight+i]*m_fHUV_Real[j*lHeight+i]; fFImag[j*lHeight+i]=fFImag[j*lHeight+i]*m_fHUV_Real[j*lHeight+i]; } //3,逆向FFT //4,得到实部 //5,逆向中心化 m_fft.ComputeConverseFFT( lWidth, lHeight, fFReal, fFImag, fTReal, fTImag, false); //把实部写到图像中 // 行 for(i = 0; i < lHeight; i++) { // 列 for(j = 0; j < lWidth; j++) { // 指向DIB第i行，第j个象素的指针 lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j; fTemp=fabs(fTReal[j*lHeight+i]); if(fTemp>255) fTemp=255; *lpSrc=(unsigned char)fTemp; } } if(fFReal) { delete []fFReal; fFReal=NULL; } if(fFImag) { delete []fFImag; fFImag=NULL; } if(fTReal) { delete []fTReal; fTReal=NULL; } if(fTImag) { delete []fTImag; fTImag=NULL; } }

 

好了!