matlab 自动生成陷波滤波器算法实现

自动扫描生成陷波滤波器

Function Code

function [aOut1,aOut2,aOut3] = optimumNotchFilter(aIn,time,D0,a,b,m)
%这仅仅是最佳陷波滤波器函数
%time(正整数，默认为1)--->生成陷波的对数，一对=两个
%D0（整数，默认11）--->陷波滤波器半径
%a（奇数，默认9）--->掩模版的长
%b（奇数，默认11）--->掩模版的宽
%m（正整数，默认为2）--->巴特沃斯滤波器的级数，越大越接近理想滤波器
%aOut1--->输出结果图
%aOut2--->输出陷波滤波器图
%aOut3--->输出原图像的噪声图

%输入参数接收及预定义
if nargin == 1
    time = 1;
    D0 = 11;
    a = 9;
    b = 11;
    m = 2;
elseif nargin == 2
    D0 = 11;
    a = 9;
    b = 11;
    m = 2;
elseif nargin == 3
    a = 9;
    b = 11;
    m = 2;
elseif nargin == 4
    b = 11;
    m = 2;
elseif nargin == 5
    m = 2;
end

%输入图像并进行预处理
a1 = double(aIn);

%对图像进行基2填充
[ra,ca] = size(a1);
maxL = max(ra,ca);
n = 1;
while(2^n<maxL)
   n = n+1; 
end
aFill = zeros(2^n,2^n);
[Ra,Ca] = size(aFill);
aFill(1:ra,1:ca) = a1;
% figure,imshow(aFill,[])%---------------Figure------------

%对原图像进行保边填充（以便在图像边缘的点不超出索引范围）
abFill = zeros(ra+floor(b/2)*2,ca+floor(a/2)*2);
[abRa,abRc] = size(abFill);
abFill(floor(b/2)+1:abRa-floor(b/2),floor(a/2)+1:abRc-floor(a/2)) = a1;
gxy = abFill;

%生成退化图像的频谱
[X,Y] = meshgrid(0:Ra-1,0:Ca-1);
aFDft = fft2(aFill.*(-1).^(X+Y));
% figure,imshow(log(abs(aFDft)+1),[])%---------------Figure------------

%生成陷波滤波器
[~,Huv1] = NotchFilter(aFDft,time,D0,1.5,m,0);
% figure,imshow(Huv1,[]);%---------------Figure------------

%采集退化图像的噪声
aNoiseDft = aFDft.*Huv1;
aNoiseTemp = ifft2(aNoiseDft);
aNoise = real(aNoiseTemp).*(-1).^(X+Y);
nxy = aNoise(1:ra,1:ca);
% figure,imshow(nxy,[]);%---------------Figure------------

%对噪声像进行保边填充
nabFill = zeros(ra+floor(b/2)*2,ca+floor(a/2)*2);
nabFill(floor(b/2)+1:abRa-floor(b/2),floor(a/2)+1:abRc-floor(a/2)) = nxy;
nxy = nabFill;

%生成最佳陷波值w,掩模版a = 33；b = 33；
aOutTemp = gxy;
for f1 = floor(b/2)+1:abRa-floor(b/2)%行的循环，看b
    for f2 = floor(a/2)+1:abRc-floor(a/2)%列的循环，看a
        Sgxy = gxy(f1-floor(b/2):f1+floor(b/2),f2-floor(a/2):f2+floor(a/2));
        Snxy = nxy(f1-floor(b/2):f1+floor(b/2),f2-floor(a/2):f2+floor(a/2));
        SgAve = mean(Sgxy(:));
        SnAve = mean(Snxy(:));
        SgnAve = mean(Sgxy(:).*Snxy(:));
        Sn2Ave = mean(Snxy(:).*Snxy(:));
        wxy = (SgnAve-SgAve*SnAve)./(Sn2Ave-SnAve.*SnAve+eps);
        aOutTemp(f1,f2) = gxy(f1,f2)-wxy*nxy(f1,f2);
    end
end
aOut1 = aOutTemp(floor(b/2)+1:abRa-floor(b/2),floor(a/2)+1:abRc-floor(a/2));
aOut2 = Huv1;
aOut3 = nxy;
% figure,imshow(aOut,[])%---------------Figure------------
end

Demo Code

clear
close all

%预处理
aIn = imread("car.jpg");
a1 = double(aIn);

%得到输入图像频谱图
[ra,ca] = size(a1);
[X,Y] = meshgrid(0:ca-1,0:ra-1);
aDft = fft2(a1.*(-1).^(X+Y));

%自定义参数
time = 1;
D0 = 9;
k0 = 2;
m = 2;
keep = 0;

%带入函数并得出陷波滤波器（巴特沃斯）
[aOut1,aOut2] = NotchFilter(aDft,time,D0,k0,m,keep);
figure
subplot(121),imshow(aOut1,[])
subplot(122),imshow(aOut2,[])

Figure

1、有零频，零频更大，显示一对陷波

%自定义参数
time = 1;
D0 = 9;
k0 = 2;
m = 2;
keep = 1;

2、有零频，零频等大，显示二对陷波

%自定义参数
time = 2;
D0 = 9;
k0 = 1;
m = 2;
keep = 1;

3、有零频，零频更大，显示5对陷波

%自定义参数
time = 5;
D0 = 9;
k0 = 2;
m = 2;
keep = 1;

4、无零频，显示4对陷波

%自定义参数
time = 4;
D0 = 9;
k0 = 2;
m = 2;
keep = 0;