在本示例中，共設計了三個函數，分別是巴特沃斯濾波器BLPF（）、巴特沃斯陷波濾波器notchFilter_BTW（）、高斯陷波濾波器notchFilter_GAUSS（）

巴特沃斯陷波濾波器參見書上6.4.5選擇性濾波器，高斯陷波濾波器參照6.3.3節中的高斯低通濾波器。

參數說明：

rows——濾波器的行數；
cols——濾波器的列數；
D0——頻域截止半徑；
n——巴特沃斯的階數；
cvtype——濾波器的數據類型和通道數，默認為雙通道浮點數

Mat BLPF(int rows,int cols,float D0,int n=1,int cvtype=CV_32FC2)
{
    Mat filt(rows,cols,cvtype,Scalar::all(
 
0));
    int cx=cols/2,cy=rows/2;
    float D02=D0*D0;
    for(int i=0;i<rows;i++)
    {
        for(int j=0;j<cols;j++)
        {
            int u=cx-j,v=cy-i;//中心坐標
            float Duv2=u*u+v*v;//距離中心半徑的平方
            float H0=1-1/(1+pow(Duv2/D02,n));
            int u1=u-130,v1=v-130;
            float D2=u1*u1+v1*v1;
            
 
float H1=1-1/(1+pow(D2/D02,n));
            //.data返回的是uchar*型指針，所以要強制轉換成浮點數型
            float* p=(float*)(filt.data+i*filt.step[0]+j*filt.step[1]);

            for(int c=0;c<filt.channels();c++)
            {
                p[c]=H1*H0;
            }

        }
    }
    return filt;
}


Mat notchFilter_BTW(
 
int rows,int cols,std::vector<cv::Point> np,
                float* D,int n=1,int cvtype=CV_32FC2)
{
    Mat filt(rows,cols,cvtype,Scalar::all(0));
    int cx=cols/2,cy=rows/2;
    int numNotch=np.size();
    float* D2=D;
    for(int i=0;i<numNotch;i++)
    {
        D2[i]=D[i]*D[i];
    }
    int uk[numNotch],vk[numNotch];//在畫面上的實際陷波坐標點
    int u_k[numNotch],v_k[numNotch];//陷波共軛點
    float Dk[numNotch],D_k[numNotch];//陷波半徑r
    float Hk[numNotch],H_k[numNotch];

    for(int i=0;i<rows;i++)
    {
        for(int j=0;j<cols;j++)
        {
            int u=cx-j,v=cy-i;//中心坐標
            for(int s=0;s<numNotch;s++)
            {
                uk[s]=u+np[s].x,vk[s]=v+np[s].y;
                Dk[s]=uk[s]*uk[s]+vk[s]*vk[s];//距離中心半徑的平方
                Hk[s]=1-1/(1+pow(Dk[s]/D2[s],n));

                u_k[s]=u-np[s].x,v_k[s]=v-np[s].y;
                D_k[s]=u_k[s]*u_k[s]+v_k[s]*v_k[s];
                H_k[s]=1-1/(1+pow(D_k[s]/D2[s],n));
            }
            //.data返回的是uchar*型指針，所以要強制轉換成浮點數型
            float* p=(float*)(filt.data+i*filt.step[0]+j*filt.step[1]);

            for(int c=0;c<filt.channels();c++)
            {
                p[c]=Hk[0]*H_k[0];
                for(int k=1;k<numNotch;k++)
                {
                    p[c]*=Hk[k]*H_k[k];
                }
            }

        }
    }
    return filt;
}
Mat notchFilter_GAUSS(int rows,int cols,std::vector<cv::Point> np,
                float* D,int cvtype=CV_32FC2)
{
    Mat filt(rows,cols,cvtype,Scalar::all(0));
    int cx=cols/2,cy=rows/2;

//    float D02=D0*D0;
    int numNotch=np.size();
    float* D2=D;
    for(int i=0;i<numNotch;i++)
    {
        D2[i]=D[i]*D[i];
    }
    int uk[numNotch],vk[numNotch];//在畫面上的實際陷波坐標點
    int u_k[numNotch],v_k[numNotch];//陷波共軛點
    float Dk[numNotch],D_k[numNotch];//陷波半徑r
    float Hk[numNotch],H_k[numNotch];

    for(int i=0;i<rows;i++)
    {
        for(int j=0;j<cols;j++)
        {
            int u=cx-j,v=cy-i;//中心坐標
            for(int s=0;s<numNotch;s++)
            {
                uk[s]=u+np[s].x,vk[s]=v+np[s].y;
                Dk[s]=uk[s]*uk[s]+vk[s]*vk[s];//距離中心半徑的平方
                Hk[s]=1-exp(-Dk[s]/(D2[s]*2));

                u_k[s]=u-np[s].x,v_k[s]=v-np[s].y;
                D_k[s]=u_k[s]*u_k[s]+v_k[s]*v_k[s];
                H_k[s]=1-exp(-D_k[s]/(D2[s]*2));
            }
            //.data返回的是uchar*型指針，所以要強制轉換成浮點數型
            float* p=(float*)(filt.data+i*filt.step[0]+j*filt.step[1]);

            for(int c=0;c<filt.channels();c++)
            {
                p[c]=Hk[0]*H_k[0];
                for(int k=1;k<numNotch;k++)
                {
                    p[c]*=Hk[k]*H_k[k];
                }
            }

        }
    }
    return filt;
}

int main()
{
    Point np[]={Point(130,130),Point(90,130),Point(130,100),Point(90,100)};//輸入陷波坐標數組
    vector<Point> vnp(np,np+4);
    float D[4]={5,10,15,20};
//    Mat filt1=notchFilter_BTW(500,600,vnp,D,2);
    Mat filt1=notchFilter_GAUSS(500,600,vnp,D);
    Mat fc1;
    extractChannel(filt1,fc1,0);
    imshow("filter ",fc1);

    waitKey();
    return 0;
}

高斯陷波濾波器的演示實例結果如下：

下面是巴特沃斯陷波濾波器的演示結果：

Opencv實現的陷波濾波器