摘要：图像的频域滤波是图像增强的一种方法。图像增强是图像处理的方法之一，有频率域法和空间域法。频率域法把图像看成一种二维信号，对其进行二维傅里叶变换的信号增强，采用低通滤波法可以去掉图像的噪声；采用高通滤波法，则可增强边缘等高频信号，使模糊的图像变得清晰。滤波是一种重要的图像预处理手段，目前，图像的频域滤波应用非常广泛，例如，基于频域滤波的织物疵点的检测、颤振信号处理中的时频域滤波方法等。
本文将主要研究常用图像滤波方法分析、频域滤波的特点及使用场合以及其如何用MATLAB实现图像的增强处理。
关键词：图像增强、频域滤波、MATLAB
目录
摘要
Abstract
1绪论-1
1.1数字图像处理的研究现状及发展方向-1
1.2研究内容和研究方法-2
1.2.1研究内容-2
1.2.2研究方法-2
2数字图像处理-2
2.1数字图像的基本概念-2
2.2数字图像处理软件-4
2.2.1 MATLAB简介-4
2.2.2其他图像处理软件的介绍-4
3图像的频率域变换-5
3.1连续傅里叶变换-5
3.1.1一维连续傅里叶变换-5
3.1.2二维连续傅里叶变换-6
3.2离散傅里叶变换-6
3.3快速傅里叶变换(FFT)-7
3.4本章小结-7
4图像增强-8
4.1图像噪声-8
4.2频率域图像增强概念-8
4.3 空间域滤波和频率域滤波的比较-9
5频率域滤波增强-9
5.1一些基本的滤波器介绍-10
5.2频率域低通滤波器-10
5.2.1理想低通滤波器-10
5.2.2巴特沃斯(Butterworth)低通滤波器-10
5.2.3高斯低通滤波器-11
5.3频率域高通滤波器-11
5.3.1理想高通滤波器-11
5.3.2巴特沃斯型高通滤波器-11
5.3.3高斯高通滤波器-12
5.4同态滤波器的简单介绍-12
5.5本章小结-12
6频域滤波增强应用实例-13
6.1理想型低通滤波器的应用-13
6.2巴特沃斯(Butterworth)型低通滤波器的应用-13
6.3高斯型低通滤波器的应用-15
6.4理想型高通滤波器的应用-16
6.5巴特沃斯型高通滤波器的应用-16
6.6高斯型高通滤波器的应用-17
6.7频域滤波应用于舰船热尾流红外图像-18
7结论-20
致谢-20
参考文献-21
附录-21

华中科技大学《数字信号分析理论实践》第六单元 信号数字滤波的概念 学习总结记录

频域滤波


对信号进行低通滤波的一种快速方法，低频信号幅度设为1，负频率镜像过去，把虚部设为0



Fs = 2048;
dt = 1.0/Fs;
T = 1;
N = T/dt;
t = linspace(0,T,N);
x1 = sin(2*pi*50*t)+sin(2*pi*300*t)+sin(2*pi*500*t);
subplot(411);
plot(t,x1);
axis([0,0.1,-2,2]);
P = fft(x1,N);
Pyy = 2*sqrt(P.*conj(P))/N;
f = linspace(0,Fs/2,N/2);
subplot(412)
plot(f,Pyy(1:N/2));
P1(N) = 0;
for k = 1:N
	P1(k) = real(P(k))+1i*imag(P(k));
end
for k = 200:N-200
	P1(k) = 0;
end
Pyy = 2*sqrt(P1.*conj(P1))/N;
subplot(413)
plot(f,Pyy(1:N/2));
x2 = ifft(P1);
subplot(414);
plot(t,real(x2));
axis([0,0.1,-2,2]);



长数据滤波



数字时域滤波器——关键得到Z变换形式的滤波器

数字差分——简单的高通滤波器

$x'(n)\approx\frac{x(n+1)-x(n-1)}{2\Delta t}\Rightarrow H(z)=az^1-az^{-1},a = \frac{1}{2\Delta t}$

a = 1/(2*dt);
For K = 1 To N-1
	x(k) = a*x(k+1)-a*x(k-1)
Next


Fs = 5000;
dt = 1/Fs;
T = 1;
N = T/dt;
t = linspace(0,T,N);
a = 5;
f = 2;
y = a*sin(2*pi*f*t)+0.3*sin(2*pi*50*t);
subplot(211)
plot(t,y);
x(N) = 0;
for i = 2:N-1
	x(i) = (y(i+1)-y(i-1))/2;
end
subplot(212)
plot(t,x);


数据平滑——简单的低通滤波器



Fs = 5000;
dt = 1/Fs;
T = 1;
N = T/dt;
t = linspace(0,T,N);
a = 5;
f = 2;
y = a*sin(2*pi*f*t)+0.8*sin(2*pi*500*t);
subplot(211)
plot(t,y);
x = y;
for i = 1:N-5
	x(i) = (y(i)+y(i+1)+y(i+2)+y(i+3)+y(i+4))/5;
end
subplot(212)
plot(t,x);



高斯平滑滤波器——将系数改进，中间值权重最大


数字滤波器设计——求取滤波器系数

时域数字滤波就是将信号和滤波器的单位脉冲响应 (Z域滤波器系数) 进行卷积分，寻找所需的滤波器的单位脉冲响应的过程称为滤波器设计，信号卷积分过程称为滤波	




注意对滤波器起始点进行处理

1、频率滤波
图像的空间域滤波：用各种模板直接与图像进行卷积运算，实现对图像的处理，这种方法直接对图像空间操作，操作简单。图像处理不仅可以在空间域进行还可以在频率域进行，把空间域的图像开窗卷积形式，变换得到频率域的矩阵点乘形式得到比较好的效果。图像频域滤波，先把图像转换到频域空间，然后对不同的频率点进行滤波，使用信号处理的技术，对图像实现滤波。比如实现图像的轮廓提取，在空间域滤波中我们使用一个拉普拉斯模板就可以提取，而在频域内，我们使用一个高通滤波模板，可以实现轮廓的提取。
图像特征与像素点数值的关系：
图像尺寸418*564，如果把每一行所有像素，一行564个点的灰度作为一维向量画图，取前三行画成三条曲线，就得到了下面的图形。
subplot(2,1,1);imshow(img), title('原始图像');line1 = img(1, :);line2 = img(2, :);line3 = img(3, :);subplot(2,1,2);hold onplot(line1, 'r');plot(line2, 'g');plot(line3, 'b');
输出结果如下图。
图像前几行的特征如下图(为了方便观察提取了不止三行)：
从图像矩阵前几行的像素特征可见，图像平坦的地方像素曲线也平坦，图像亮的地方就是图像像素剧烈变化的地方，图像像素值发生较大差异的地方也是图像发生突变的地方，这些位置一般就是图像轮廓。图像的频率体现了图像中灰度变化剧烈程度，是灰度在平面空间上的梯度。所以前面讲的空域滤波，可以使用平滑来滤除噪声实现平滑，从上面曲线图像上看，可以按照信号处理思想来理解，平滑滤波就是频率低通。对应的，频率高通滤波就是空域的提取边界。
2、图像傅里叶变换
傅里叶变换可以将一个时域信号转换成在不同频率下对应的振幅及相位，其频谱就是时域信号在频域下的表现，而反傅里叶变换可以将频谱再转换回时域的信号。
imgPath = 'E:opencv_picsrc_picpic4.bmp';img = imread(imgPath);img=rgb2gray(img); f=fft2(im2double(img)); %FFTF=f; %FFT频谱T=log(F+1); %频谱对数变换subplot(1,2,1),imshow(img),title('原始图像');subplot(1,2,2),imshow(T,[]),title('原始图像其频谱图');
输出结果如下，幅度图在四个角出现亮光。因为实信号以fs为采样速率的信号在 fs/2处混叠，所以实信号fft的结果中前半部分对应[0, fs/2]，后半部分对应[ -fs/2, 0]。横向和纵向都需要把频率转换到[ -fs/2, fs/2]的区间，可以使用fftshift函数。
在数字图像处理中，常常需要将F(u,v)的原点移到N*N频域的中心，以便能清楚地分析傅里叶谱的情况，平移前空域、频域原点均在左上方。而fftshift的作用就是这样，将0频谱移到正中心。
Y = fftshift(X) 通过将零频分量移动到数组中心，重新排列傅里叶变换 X。
如果 X 是向量，则 fftshift 会将 X 的左右两半部分进行交换。
如果 X 是矩阵，则 fftshift 会将 X 的第一象限与第三象限交换，将第二象限与第四象限交换。
如果 X 是多维数组，则 fftshift 会沿每个维度交换 X 的半空间。
a =     1     2     3     4     6     7     8     9>> fftshift(a)ans =     8     9     6     7     3     4     1     2
平移后的图：
傅里叶变换的能量集中在频率很小的圆内，当D0增大时能量衰减很快，高频部分虽然携带的能量很少，但是包含丰富的边界和细节信息，所以当截止频率D0变小时，虽然亮度足够(因能量损失不大)，但图像变模糊。
3、频域滤波
傅里叶变换可以把图像从空域变换到频域，而傅里叶反变换可以将图像的频谱变换为空域图像。可以利用图像空域和频域之间的对应关系，尝试将空域卷积滤波变换为频域滤波，而后再将频域滤波处理后的图像反变换回空间域，从而达到图像增强的目的。
G=imnoise(img,'gaussian', 0, 0.05);%模拟均值为0方差为0.05的高斯噪声，H=fft2(im2double(G)); %FFTH=fftshift(H); %FFT频谱平移T=log(abs(H)); %频谱对数变换subplot(2,2,1), imshow(G),title('添加高斯噪声图像');subplot(2,2,2),imshow(T, []),title('频谱图');
使用傅里叶变换得到图像的频域幅度，应用傅里叶反变换得到空域图像.
img2 = (ifft2(ifftshift(H))); %===频域的图反变换到空域img3 = im2uint8(mat2gray(img2)); %===取其灰度图subplot(2,2,3),imshow(img3);title('anti-Fourier');
反变换效果如下图。
那么就可以在频率对图像进行滤波，然后反变换回空域就得到滤波后的图像。频率图可见在图像中心，是频率为0位置，所以如果限定只有低频通过则为低通，如果抑制频谱图中心位置，则是高通。可见图像在频域滤波是非常直观的。
1)低通滤波器：
其中D(u,v)为频率域上(u,v)点到中心的距离，D0由自己设置.
2)高通滤波器:
高通滤波器同低通滤波器非常类似，只不过二者通过的波正好是相反的.
3)高斯低通滤波
4)高斯高通滤波
高斯频率低通滤波：
S=fftshift(fft2(G)); %G是添加高斯噪声的图像[M,N]=size(S);                                          d0=5; %GLPF滤波，d0=5，15，30                    n1=floor(M/2);                          n2=floor(N/2);                           for i=1:M     for j=1:N        d=sqrt((i-n1)^2+(j-n2)^2);                 h=1*exp(-1/2*(d^2/d0^2));          S(i,j)=h*S(i,j);                       endend S=ifftshift(S);                           S=uint8(ifft2(S));    subplot(2,2,4),imshow(S),title('高斯低通滤波图像');
输出图像如下：
高斯滤波函数标准差变大，则高斯函数的傅里叶变换图像频率范围变宽，允许通过的频率范围变大，对应空域允许通过的图像变换更加多，可以允许的噪声更多。相反，如果标准差变小，频域范围变窄，允许通过的频率范围变小，对应空域对低频要求更加严格被抑制的高频更多，所以图像更加平滑，可能出现模糊现象。下面展示了不同的σ标准差情况下的滤波效果。
改为高通滤波，结果如下图：

使用Matlab对给定的信号和预期需求，设计了几种连续时间滤波器，包括四种最常见的连续滤波器：巴特沃斯(Butterworth)、切比雪夫滤波器I型和II型(chebyshev)、考尔滤波器(Cauer)；验证了滤波器性能，并进行了分析和对比，给出了代码，并将滤波器设计过程常用操作简单封装成函数，简化编程。
待处理信号
这里使用gamme.m函数，生成一段依次从Do、 Re、Mi、...一直到升Do的音频，函数输入是每个频率保持的长度 和 采样频率，输出是一段二维数组，即生成的音频。
function [gam,t] = gamme(duree,fe)

% [do re mi fa sol la si doo]
freqnotes = [262 294 330 349 392 440 494 523];
t = 0:1/fe:duree;

gam = [];
for note=1:8, 
    gam = [gam, sin(2*pi*freqnotes(note)*t)];
end
N = length(gam);
t = (0:N-1)/fe;
调用函数，生成音频，并使用matlab的soudnsc函数试听
clear all
fe = 8192 ; %采样频率
duree = 1; %每个音符持续1秒
[sig,t] = gamme(duree,fe);
soundsc(sig, fe); % listen
对原始信号进行频率分析
fmin = 250;
fmax = 550;

f = [0:length(sig)-1]*(fe/length(sig));
S = fft(sig)/fe;

figure()
plot(f, abs(S));
axis([fmin fmax 0 0.5])
连续时间滤波器
滤波器指标
我们希望除去后三个音，即La So 升Do，对照这个表，得到滤波器指标
类型：低通
通带容许的最小增益：-3dB
阻带增益：-40dB
过渡带中心频率：420Hz （在Sol和La之间）
过渡带宽度：100Hz （这里100Hz有些宽，为的是让不同的滤波器的对比更加明显，我们会在后面看到它的影响）
之所以通带并不是完全通过，是因为大多滤波器会有抖动存在，并且即使没有抖动，滤波器的频域响应也不会是水平直线，因此设计时应当留有一定的裕度
滤波器设计
1、常量定义
Ap = 3;
Aa = 40; %这里的增益不用加负号
Deltaf = 100;
fc = 420;
fa = fc + Deltaf/2;
fp = fc - Deltaf/2;
2、计算滤波器阶数
[nc1 wnc1]= buttord(2*pi*fp,2*pi*fa,Ap,Aa,'s')
[nc2 wnc2]= cheb1ord(2*pi*fp,2*pi*fa,Ap,Aa,'s')
[nc3 wnc3]= cheb2ord(2*pi*fp,2*pi*fa,Ap,Aa,'s')
[nc4 wnc4]= ellipord(2*pi*fp,2*pi*fa,Ap,Aa,'s') %cauer滤波器
- 结果分析
nc代表需求的阶数，wnc代表对应滤波器的【由通带到过度带】的频率，这个频率对于不同的滤波器来说是不一样的，直观上可以理解为，这是因为每个滤波器的频域响应的”形状“都不一样，因此它们与【-3dB的水平线】的交点的横坐标不一样。
命名的‘c’代表连续时间，即模拟滤波器
此外，可以看到，Butterworth的阶数最高，是20，两个Chebyshev阶数相同，是8，而Cauer滤波器只需要5阶就可以
通常情况下，滤波器的阶数对应了滤波器的复杂程度和成本，阶数越高，成本越高，为了达到相同的性能指标，Butterworth需要更高的阶数，Cauer需要最少的，但是有失就有得，Butterworth的好处在于它在通带和阻带都没有抖动，是光滑的，Chebyshev I型在通带有抖动，Chebyshev II型在阻带有抖动，Cauer在通带和阻带都有抖动。
这体现了一种tradeoff的思想，如果想要滤波器在通带和阻带没有抖动，那么代价就是需要更高的滤波器阶数，如果想要更简单的滤波器结构，即阶数低，那么代价是滤波器的频域响应会有抖动
这个抖动的影响，我们会在后面的仿真中看到
3、计算低通滤波器参数
[numc1, dend1] = butter(nc1,wnc1,'low','s');
[numc2, dend2] = cheby1(nc2,Ap,wnc2,'low','s');
[numc3, dend3] = cheby2(nc3,Aa,wnc3,'low','s');
[numc4, dend4] = ellip(nc4,Ap,Aa,wnc4,'s');
值得注意的是，四个函数的参数是不同的，butter默认只需要输入【阶数，ord函数计算出来的频率】，但cheby1需要输入【阶数，通带增益，ord函数计算出来的频率】，cheby2则是把通带增益换成了阻带增益，ellip，即cauer滤波器则需要两个都输入。这与滤波器的频域响应抖动是对应的：Cheby1 通带有抖动，Cheby2阻带有抖动，Cauer 通带阻带都有抖动。
对比四种滤波器的频域响应
f = linspace(fmin*0.1,fmax*2,1000);
Hc1 = freqs(numc1,dend1,2*pi*f);
Hc2 = freqs(numc2,dend2,2*pi*f);
Hc3 = freqs(numc3,dend3,2*pi*f);
Hc4 = freqs(numc4,dend4,2*pi*f);
figure()
hold on
plot(f,20*log10(abs(Hc1)))
plot(f,20*log10(abs(Hc2)))
plot(f,20*log10(abs(Hc3)))
plot(f,20*log10(abs(Hc4)))
axis([fmin*0.1 fmax*2 -70 5])
grid on
legend('Butterworth',' Chebyshev de type 1' ,' Chebyshev de type 2' ,' Cauer');
除了抖动，还可以看到过渡带（从通带到阻带的部分）的宽度也不同，紫色和橙色下降地比较迅速，过渡带比较窄，蓝色和黄色下降地比较缓慢，过渡带比较宽，理论上来讲，过渡带越窄，滤波器越接近理想的滤波器模型，滤波性能越好。
应用滤波器处理信号
使用lsim命令得到滤波后的信号
sigf1 = lsim(tf(numc1,dend1),sig,t);
sigf2 = lsim(tf(numc2,dend2),sig,t);
sigf3 = lsim(tf(numc3,dend3),sig,t);
sigf4 = lsim(tf(numc4,dend4),sig,t);
性能分析
可以看一下得到的信号的频域与原信号有什么不同
频域上
figure()
hold on
subplot(3,2,1)
FFTplot(sig,fe,fmin,fmax) %这里的FFT是手写函数，用于绘制信号的频谱，很简单，详细定义见文末
subplot(3,2,3)
FFTplot(sigf1,fe,fmin,fmax)
subplot(3,2,4)
FFTplot(sigf2,fe,fmin,fmax)
subplot(3,2,5)
FFTplot(sigf3,fe,fmin,fmax)
subplot(3,2,6)
FFTplot(sigf4,fe,fmin,fmax)
时域上
figure()
hold on
subplot(3,2,1)
plot(t,abs(sig))
title('Raw signal')
subplot(3,2,3)
plot(t,abs(sigf1))
title('Butterworth')
subplot(3,2,4)
plot(t,abs(sigf2))
title('Chebyshev de type 1')
subplot(3,2,5)
plot(t,abs(sigf3))
title('Chebyshev de type 2')
subplot(3,2,6)
plot(t,abs(sigf4))
title('Cauer')
原始信号每秒依次发出Do，Re，Mi....幅值都是1，所以是个大方块
频域和时域的图很像，这是因为这个例子的特殊性：时间上多个频率是独立的，即这段信号中，每个频率的成分都是依次出现的，不存在多个频率同时发生的情况。
可以看到：
第七个和第八个音都基本被很好的消除了，第六个音差一点
我们设计阻带的增益是-40dB，根据Bode的计算公式，这意味着“滤波后的该频率的幅值是原先的 
 倍”
第五个音并不是我们想要消除的，但是都受到了不同程度的影响，这是因为我们的滤波器指标有误，其中的过渡带太宽了，导致第五个音受到了影响
不同滤波器在通带的抖动情况也不一样，Butter和Cheby2在通带无抖（但依旧在第四个音有可见的衰减！），Cheby1和Cauer在通带有抖，这影响了最终得到的信号质量
重新设计性能指标
从上面的仿真可以看到，因为指标的错误，第五个音被意外地滤掉了一大部分，并且四个滤波器的通带表现也不好，要么是有抖，要么是有可见的衰减，为了得到更好的滤波效果，我们更改滤波器指标，让它更加严格。
Ap = 1; %通带增益从-3到-1
Aa = 40;
Deltaf = 20; %过渡带带宽从100到20
fc = 420;
fa = fc + Deltaf/2;
fp = fc - Deltaf/2;
重新计算四种滤波器，可以得到性能更加好但阶数也更高的滤波器
值得一提的是，此时Butterworth滤波器的系数计算结果会发散，我怀疑这意味Butterworth不能满足对应的频域指标
我也把计算四种连续时间低通滤波器并画图的代码封装成了函数，只需要一行，就可以得到下图，函数放在了文末
[numc1,denc1]=FiltrecoePBA(Ap,Aa,fc,Deltaf,'c1',1); %c1代表cheby1型，最后的1代表输出频域对比图
新得到的这几种滤波器都能更好的完成滤波任务，这体现在：第五个音没有被滤掉，后三个音都被很好地滤掉。
手动封装的函数
为了便于对比各个滤波器的性能，让程序更加简洁，封装了几个常用操作：
输出信号傅里叶图谱
输入预期指标，计算四种连续低通滤波器，并输出四种滤波器频域对比图
输入信号和滤波器，输出滤波后信号，并显示滤波前后频域分析
信号频域分析
函数说明和例程
Draw the Fourier transform of a signal.
Inputs
sig: signal
fe: signal sampling frequency
Optional input
[min_ max _] represents the frequency domain interval.
-------------------------------------------------------------------------
example: %复制之后直接跑即可
    fe = 1000;            % Sampling frequency                    
    T = 1/fe;             % Sampling period       
    L = 7000;             % Length of signal
    t = (0:L-1)*T;        % Time vector
    sig = 0.7*sin(2*pi*50*t) + sin(2*pi*128*t) + 0.2*randn(size(t)); 
    FFTplot(sig,fe)
例程中生成了一段带有噪声的两个正弦信号的叠加，输出是该段信号的频域图
函数正文
function no_output = FFTplot(sig, fe,min_, max_)

if nargin == 2
    min_ = 0 ;
    max_ = fe/2;
elseif nargin < 2
    disp('input error!')
    return 
end

f = [0:length(sig)-1]*(fe/length(sig));
S = fft(sig)/fe;
plot(f, abs(S),'LineWidth',1.5);

axis([min_ max_ 0 max(abs(S))*1.2])
grid on

end
计算四种连续低通滤波器
函数说明和例程
Calculating Transfer Function Parameters for Four Continuous-Time Low-Pass Filters
Input: 
        Ap,  ondulation bande passante  
        Aa,  attenuation bande attenuee
        fc, 
        Deltaf,
        type, 'b'/'c1'/'c2'/'e'
        plots = 1 : show bode plot  
                else : without plot 
output: num/den of filter
-----------------------------------------------------------------------------------------------
example: %复制之后直接跑即可
        Ap = 3;
        Aa = 40;
        Deltaf = 100;
        fc = 420;
        type = 'c1'
        plots = 1
        [numc,denc] = FiltrecoePBA(Ap,Aa,fc,Deltaf,type,plots); % calculate filter
例程定义了对应的低通滤波器性能指标，调用该函数计算cheby1类型滤波器，返回滤波器分子分母参数，并画出四种滤波器频域响应
函数正文
function [numc,denc] = FiltrecoePBA(Ap,Aa,fc,Deltaf,type,plots)
fa = fc + Deltaf/2;
fp = fc - Deltaf/2;

[nc1 wnc1]= buttord(2*pi*fp,2*pi*fa,Ap,Aa,'s');
[nc2 wnc2]= cheb1ord(2*pi*fp,2*pi*fa,Ap,Aa,'s');
[nc3 wnc3]= cheb2ord(2*pi*fp,2*pi*fa,Ap,Aa,'s');
[nc4 wnc4]= ellipord(2*pi*fp,2*pi*fa,Ap,Aa,'s');

[numc1, denc1] = butter(nc1,wnc1,'low','s');
[numc2, denc2] = cheby1(nc2,Ap,wnc2,'low','s');
[numc3, denc3] = cheby2(nc3,Aa,wnc3,'low','s');
[numc4, denc4] = ellip(nc4,Ap,Aa,wnc4,'s');



if type == 'b'
    numc_ = numc1;
    denc_ = denc1;
elseif type == 'c1'
    numc_ = numc2;
    denc_ = denc2;
elseif type == 'c2'
    numc_ = numc3;
    denc_ = denc3;
elseif type == 'e'
    numc_ = numc4;
    denc_ = denc4;
else 
    disp(type)
    disp('not define!')
    return 
end

if sum(isnan(numc_))>0 || sum(isnan(denc_))>0
    disp(['Error: filter of type ',type,' diverge!'])
    return 
else 
    numc = numc_;
    denc = denc_;
end




if plots==1
    f = linspace(fc/10,fc*4,1000);
    Hc1 = freqs(numc1,denc1,2*pi*f);
    Hc2 = freqs(numc2,denc2,2*pi*f);
    Hc3 = freqs(numc3,denc3,2*pi*f);
    Hc4 = freqs(numc4,denc4,2*pi*f);
    figure()
    hold on
    plot(f,20*log10(abs(Hc1)))
    plot(f,20*log10(abs(Hc2)))
    plot(f,20*log10(abs(Hc3)))
    plot(f,20*log10(abs(Hc4)))
    axis([fc/10,fc*4 -2*Aa 5])
    grid on 
    legend('Butterworth',' Chabyshev type 1' ,' Chabyshev type 2' ,' Cauer');
    hold off
end



end
输入信号和滤波器，输出滤波后信号，并显示滤波前后频域分析
函数说明及例程
Given a signal and a continuous filter, draw the spectrum of the raw signal 
and the filtered spectrum, return the filtered signal
Inputs
sig: Signal
fe: Sampling Frequency
numc : the molecular coefficient of the continuous filter transfer function.
denc: The denominator of the continuous filter transfer function.
Optional input
min_ max_ : Frequency Interval
Ouput = filtered signal
----------------------------------------------------------------------------------------
Example:  %复制之后直接跑即可
    fe = 1000;            % Sampling frequency                    
    T = 1/fe;              % Sampling period       
    L = 7000;             % Length of signal
    t = (0:L-1)*T;        % Time vector
    sig = 0.7*sin(2*pi*50*t) + sin(2*pi*128*t) + 0.2*randn(size(t));  % generate signal 
    Ap = 3; Aa = 40; Deltaf = 50; fc = 128; type = 'c1';, plots = 0;
    [numc,denc] = FiltrecoePBA(Ap,Aa,fc,Deltaf,type,plots); % calculate filter
    FiltreA(sig,fe,numc,denc,30,150);
例程中生成了一段带有噪声的两个正弦信号的叠加，并生成了一个c1类型的低通滤波器，输出滤波后的信号，并显示滤波前后的频域对比
函数正文
function [sigf] = FiltreA(sig, fe, numc, denc, min_, max_)

t = [0:1/fe:1/fe*(length(sig)-1)];   
sigf = lsim(tf(numc,denc),sig,t);

subplot(2,1,2)
FFTplot(sigf, fe,min_, max_)
title('Signal filtered')
grid on 
subplot(2,1,1)
FFTplot(sig, fe,min_, max_)
grid on 
title('Raw Signal')

end