频率域滤波去除周期性噪声

傅里叶变换后的频率域去噪

（做些小小更改，让变换结果更加清晰合理）（2021年1月1日17:36:36）
去除周期性波纹噪声最重要在于
1.频率域变换问题关键在于如何准确找到噪声点的位置。这里可以用类似矩阵扫描的方法找出某个点，其满足大于其上下左右各点的值（找到局部极大值点），同时满足大于某个阈值，我给定的是大于图像均值（中心点亮度）的4/5左右，即可确定准确的坐标位置。进而用巴特沃斯滤波进行处理。
2.确定合适的参数：截止半径r(相对于噪声点中心）。即自定义函数function d = findd0(coor,image)。大概是噪声中心到频谱图中心的距离，阈值：T，截止半径：R（相对于原点）
3.运用合适的滤波，对于不同的噪声用不同的滤波公式会有不同的效果，选择合适的滤波会有更好的效果。

功能函数代码部分：(matlab)

%频率域变换去除周期性噪声
function re = dpressnoise(img)
t1=myfft(img);%一维横向变换
f=myfft(t1.');%纵向变换
[m,n]=size(f);%获取傅里叶变换后频谱图的大小
[M,N]=size(img);%获取原来图像大小
f=myfftshift(f);      %使图像对称,中心化
% r=real(f);          %图像频域实部
% i=imag(f);          %图像频域虚部
Margin=log(sqrt((real(f)).^2+(imag(f)).^2));      %图像幅度谱，加log便于显示
phase=real(angle(f)*180/pi);     %图像相位谱
figure(1)
subplot(1,3,1),imshow(linear1(img)),title('源图像');%拉伸显示
subplot(1,3,2),imshow(Margin,[]),title('图像幅度谱');
subplot(1,3,3),imshow(phase,[]),title('图像相位谱');
T=5/6*max(max(Margin));%选取频谱图最亮点的5/6作为阈值，可以根据实际情况进行修改
R=10;%频谱图中心亮点及其周围亮点是影像图像亮度及信息不是噪声，设置一个截止半径，使得中心数据不被破坏
coor=[];%设置一个变量存储检测到的噪声点的坐标信息
for i=1:m-2
    for j=1:n-2
        if(Margin(i+1,j+1)>T&&sqrt((m/2-i)^2+(j-n/2)^2)>R&&Margin(i+1,j+1)>Margin(i+1,j)&&Margin(i+1,j+1)>Margin(i+1,j+2)&&Margin(i+1,j+1)>Margin(i,j+1)&&Margin(i+1,j+1)>Margin(i+2,j+1))
            coor=[coor;i+1,j+1];
        end
    end
end
d0=findd0(coor,Margin);%计算每个噪声点对应到中心原点的距离
H=createH(d0,coor,Margin);%计算整个频谱图对应的滤波矩阵
f=f.*H;%将滤波矩阵与频谱图相乘
%逆变换
f=myfftshift(f);
g=myifft(f);
k=myifft(g.');
k=k(1:M,1:N);%截取原图像大小的图像，由于傅里叶变换中补零导致会出现黑色补零区域，现在剪裁掉
figure(2)
imshow(linear2(k)),%拉伸显示原图像
title('去噪处理结果')
re=k;
end

function H=createH(d0,coor,image)%创建H矩阵，coor指给定的噪声点中心坐标，image指频谱图像
count=size(coor,1);
[m,n]=size(image);
d=zeros(m,n,count/2);%申请矩阵
H=ones(m,n);
D0=zeros(count/2,1);
for i=1:m
    for j=1:n
        for k=1:count/2
        d(i,j,k)=sqrt((i-coor(k,1))^2+(j-coor(k,2))^2)*sqrt((i-coor(count-k+1,1))^2+(j-coor(count-k+1,2))^2);
        D0(k,1)=d0(k,1)^2;
        end
        for k=1:count/2
            if(d(i,j,k)==0)
            H(i,j)=0;
            else
            H(i,j)=H(i,j)/(1+(D0(k,1)/d(i,j,k))^2);
            end
        end
    end
end
end

function d = findd0(coor,image)%根据给出的坐标大概确定截止半径,coor指给定的噪声点中心坐标
[M,N]=size(image);
k=size(coor,1);
dis=[];
for i=1:k%循环坐标个数
    x=coor(i,1);%获取x坐标
    y=coor(i,2);%获取y坐标
    while(image(x,y+1)<image(x,y))%找出距最近的极小点
        y=y+1;
    end
    dis=[dis;M/2+1-coor(i,1)];%算出距离保存
end
d=dis;
end

myfft代码：

#%快速傅里叶变换，不够2的整数幂的个数，末尾自动补齐0
function ret_val = myfft(Vector)
%因为输入的数据可能不是2的整数次幂，变换使得计算更加方便
[m,n]=size(Vector);%输入信号矩阵大小
num=ceil(log2(n));%向上取整
N=2^num;
vector=zeros(m,N);%申请足够大小矩阵
vector(:,1:n)=Vector(:,:);%将变换后的信号输入，不足补零
%变址运算
% temp=zeros(1,N);
% for p=1:m
%     for q=0:N-1
%         t=bin2dec(fliplr(dec2bin(q,num)));%变址运算
%         temp(1,q+1)=vector(p,t+1);
%     end
%     vector(p,:)=temp(1,:);
% end
%变址运算
for line=1:m%循环行数
    j1 = 0;
    for i = 1 : N%循环个数
        if i < j1 + 1
            tmp = vector(line,j1 + 1);
            vector(line,j1 + 1) = vector(line,i);
            vector(line,i) =tmp;
        end
        k = N / 2;
        while k <= j1
            j1 = j1 - k;
            k = k / 2;
        end
        j1 = j1 + k;
    end
end
for lines=1:m%循环行数
    for L=0:num-1%L表示运算等级或者层数
        dis=2^L;%dis表示奇偶组之间的距离
        for id=1:2^(num-L-1) %循环当前层数组数
            %进行同址运算
           for idx=1:dis%循环组内个数
               x1=(id-1)*2*dis+idx;%求得奇数数组的索引值
               x2=(id-1)*2*dis+dis+idx;%对应偶数数组的索引值
            temp1=vector(lines,x1)+vector(lines,x2)*W(L,(x1-1));%中间变量保存相应奇偶数组数据
            temp2=vector(lines,x1)-vector(lines,x2)*W(L,(x1-1));
            vector(lines,x1)=temp1;%所谓存入之前地址
            vector(lines,x2)=temp2;
           end
        end
    end
end
ret_val =vector;
function val=W(L,x)%旋转因子当层数为L,索引值为x
val=exp(-1j*2*pi*x/2^(L+1));
end
end

myfftshift代码

%我的fftshift函数
function re=myfftshift(a)
[m,n]=size(a);
m1=ceil(m/2);%向上取整
n1=ceil(n/2);
%左右半边进行交换
temp1=zeros(m,n1);%申请矩阵空间
temp2=zeros(m,n-n1);
temp1(:,:)=a(:,1:n1);%存储左半数据
temp2(:,:)=a(:,n1+1:n);%存储右半数据
a(:,1:n-n1)=temp2;%赋值给左半边
a(:,n-n1+1:n)=temp1;%赋值给右半边
%上下半边进行交换
temp1=zeros(m1,n);%申请矩阵空间
temp2=zeros(m-m1,n);
temp1(:,:)=a(1:m1,:);%存储上半数据
temp2(:,:)=a(m1+1:m,:);%存储下半数据
a(1:m-m1,:)=temp2;%赋值给上半边
a(m-m1+1:m,:)=temp1;%赋值给下半边
re=a;
end

运行截图：
原图像

幅度谱图：

相位谱图：

去噪后图像：

参考文献：
论文：图像线状和网格状噪声的去除方法
陶胜