1 理想低通滤波

%理想低通
I = imread('fig.png');
I=rgb2gray(I);
figure(1);
subplot(221),imshow(I);
title('原图像');
I=imnoise(I,'gaussian');%%加入高斯白噪声
subplot(222),imshow(I);
title('加入噪声后的图像');
s=fftshift(fft2(I));
subplot(223), imshow(log(abs(s)),[]); 
title('图像傅里叶变换取对数所得频谱');
[a,b]=size(s);
a0=round(a/2);
b0=round(b/2);
d=50;
for i=1:a 
    for j=1:b 
        distance=sqrt((i-a0)^2+(j-b0)^2);
        if distance<=d
            h=1;
        else
            h=0;
        end
        s(i,j)=h*s(i,j);
    end
end
s=uint8(real(ifft2(ifftshift(s))));
subplot(224),imshow(s);
title('理想低通滤波所得图像'); 

2 高斯低通滤波

%高斯低通
I=rgb2gray(imread('fig.png'));
IA=imnoise(I,'gaussian');%%加入高斯白噪声
subplot(222),imshow(IA);
title('加入噪声后的图像');
[f1,f2]=freqspace(size(IA),'meshgrid');
D=100/size(IA,1);
r=f1.^2+f2.^2;
Hd=ones(size(IA));
for i=1:size(IA,1)
    for j=1:size(IA,2)
        t=r(i,j)/(D*D);
        Hd(i,j)=exp(-t);
    end
end
Y=fft2(double(IA));
Y=fftshift(Y);
Ya=Y.*Hd;
Ya=ifftshift(Ya);
Ia=real(ifft2(Ya));
subplot(221),imshow(uint8(I));
title('原图像');
subplot(223),imshow(uint8(Ia));
title('高斯低通滤波处理');

3 巴特沃斯低通滤波

clear all;
%巴特沃斯低通
figure(2);
I=rgb2gray(imread('fig.png'));
subplot(221),imshow(I);
title('原始图像');
I1=imnoise(I,'gaussian');%%加入高斯白噪声
subplot(222),imshow(I1);
title('加入噪声后的图像');
f=double(I1);%数据类型转换
g=fft2(f);%图像傅里叶转换 
g=fftshift(g);%傅里叶变换平移
F2=log(abs(g));%对傅里叶变换结果取绝对值，然后取对数
[N1,N2]=size(g);%傅里叶变换图像尺寸
n=2;%参数赋初始值
d0=30;
n1=fix(N1/2);%数据圆整
n2=fix(N2/2);%数据圆整
for i=1:N1%遍历图像像素
    for j=1:N2 
        d=sqrt((i-n1)^2+(j-n2)^2);
        if d==0 
            h=0; 
        else
            h=1/(1+(d/d0)^(2*n));
        end
        result(i,j)=h*g(i,j);%?图像矩阵计算处理
    end
end
F3=log(abs(result));%对傅里叶变换结果取绝对值，然后取对数
result=ifftshift(result);
X2=ifft2(result);
X3=uint8(real(X2));
subplot(223),imshow(uint8(X3))
title('Butterworth低通滤波图像');

此处是MATLAB 的代码

wp=2*pi*5000;ws=2*pi*12000;Rp=2;Rs=30;
[N,wc]=buttord(wp,ws,Rp,Rs,'s');% s is the anology filter 
[B,A]=butter(N,wc,'s');
sk=0:511; fk=0:14000/512:14000;wk=2*pi*fk;
Hk=freqs(B,A,wk);
plot(fk/1000,20*log10(abs(Hk)),'r') ;hold on
xlabel('频率(khz)');ylabel('幅度（db）');
axis([0,14,-40,5]);

[N,wc]=buttord(wp,ws,10,Rs,'s');% s is the anology filter 
[B,A]=butter(N,wc,'s');
Hk=freqs(B,A,wk);
plot(fk/1000,20*log10(abs(Hk)),'g');

[N,wc]=buttord(wp,ws,20,Rs,'s');% s is the anology filter 
[B,A]=butter(N,wc,'s');
Hk=freqs(B,A,wk);
plot(fk/1000,20*log10(abs(Hk)),'b');

legend('Rp=2','Rp=10','Rp=20','Location','northeast','Orientation','horizontal')
 

 运行结果： 分析，Rp 越大，通达的振幅越大

 参考“：

用matlab做巴特沃斯低通滤波器

趁着暑假，做一个心电图的matlab实验，遇到了滤波器问题，网上代码比较杂乱，做了一个汇总整理。
主要做了一个简单的低通滤波器并以三角函数为例子进行低通滤波。

基本数据

fk=100  %采样频率
N=1024  %采样个数
n=-N/2:N/2-1
f=n*fk/N
t=n/fk
y=sin(2*pi*10*t)

当然 f也可以表示成 f=linspace(-fk/2,fk/2,N);

第一步：做出原信号的频谱函数

注意：用fft函数作频谱分析，得到的是0~fk内的频谱
而用fftshift函数得到-fk/2~fk/2内的频谱

fft_y=fft(y);
fftshift_y=fftshift(fft_y);
f=linspace(-50,50,1024);
plot(f,abs(fftshift_y));

第二步：做出巴特沃斯低通滤波器

用buttord函数求出阶数以及Wn

butter函数是求Butterworth数字滤波器的系数，在求出系数后对信号进行滤波时用filter函数
[g,Wn]=buttord(Wp,ws,Rp,Rs);
Wp=fp/(0.5fk)=20/50=0.4(通带截至频率)
Ws=fs/(0.5*fk)=30/50=0.6(阻带截止频率)
Rp=1; 通带最大衰减
Rs=30; 阻带最小衰减
代码如下：

[g,Wn]=buttord(0.4,0.6,1,30);

用butter 求出差分方程的系数b a

不懂差分方程及filter用法的朋友可以点以下链接，参考这位大佬的理解
差分方程介绍

[b,a]=butter(g,Wn);

至此 准备工作都做完了 ，来看看这滤波器的相频特性啥样子

滤波器的相频特性

 [q,w]=freqz(b,a,256);
 plot(w*fs/(2*pi),abs(q))

所以说 咱写出来的滤波器长这个样子
20Hz为通带截至频率
30Hz为阻带截至频率

第三步：用滤波器过滤信号 并得出频谱图

用filter过滤

用filter函数就能把原信号y过滤成k函数 ，需要注意的是，它们都是时域函数，需要进行fft变换才能看到频谱图

 k=filter(b,a,y);
 fft_k=fft(k);
 fftshift_k=fftshift(fft_k);
 plot(f,abs(fftshift_k));

过滤后的结果如下：

没错，没变！我们低通滤波器允许20Hz以下的信号不变，而我们原信号的频率就是10Hz！

对高频的信号的低通滤波

这里就不赘述了，对y=sin(2pi40*t)的信号，依然用上面的低通过滤的结果如下图：

不要看它长的怪，注意一下幅值，约等于0，即把40Hz的信号过滤掉了！
行文至此，点个赞呗亲！