MATLAB滤波器设计

MATLAB有专门的滤波器设计与分析的工具箱fdatool（Filter Design and Analysis Tool），里面可以很方便地设计各种各样的滤波器，这里介绍下如何简单设计滤波器并进行使用。
在命令栏中输入fdatool
PS：fdatool在后续的版本中会换成filterDesigner
可以打开fdatool的界面如下所示

在界面的左上角是关于设计滤波器的一些信息，包括设计滤波器的结构、阶数，稳定性和设计状态。

右上角是关于滤波器的幅频响应（默认），通过Analysis可以得到其他曲线，如相频响应、群延时、零极点等等。

界面的下方是关于滤波器的具体参数，包括响应类型、设计方法、采样率、滤波器阶数、截止频率、通带，阻带波动等等

可以看出有低通、高通、带通、带阻类型的滤波器，以及是FIR还是IIR滤波器，具体类型及数值可以根据设计的要求进行选定，这里就不再叙述。
下面给出一个简单的例子。
对于信号
$$y=s_1\text{+}{s}_2+s_3$$
其中，$s_1=\sin \left(2\pi f_{1}t\right)$，$s_2=\sin \left(2\pi f_{2}t\right)$，$s_3=\sin \left(2\pi f_{3}t\right)$，$f_1=60\text{Hz}$ ，$f_2=130\text{Hz}$，$f_3=1000\text{Hz}$，要求设计一个低通滤波器，滤除频率130Hz和1000Hz的信号。
我们可以按照要求很轻易地设计出一个滤波器

为了能够调用这个滤波器，我们可以将其保存成一个函数

如myfilter.m，这样后续就可以直接调用，但是这种调用方式相当于在运行到该命令时重新设计滤波器，对于一些阶数比较大的滤波器可能会耗时比较长。保存的函数为

function Hd = myfilter
%MYFILTER Returns a discrete-time filter object.

% MATLAB Code
% Generated by MATLAB(R) 9.1 and the Signal Processing Toolbox 7.3.
% Generated on: 07-Feb-2021 11:37:07

% Butterworth Lowpass filter designed using FDESIGN.LOWPASS.

% All frequency values are in Hz.
Fs = 10000;  % Sampling Frequency

Fpass = 80;          % Passband Frequency
Fstop = 90;          % Stopband Frequency
Apass = 0.1;         % Passband Ripple (dB)
Astop = 80;          % Stopband Attenuation (dB)
match = 'stopband';  % Band to match exactly

% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.lowpass(Fpass, Fstop, Apass, Astop, Fs);
Hd = design(h, 'butter', 'MatchExactly', match);

% [EOF]

也可以直接导出滤波器的系数，这样就可以直接调用系数进行滤波，会加快程序运行的速度。具体在：File——Export，界面如下：

可以看出，系数可以导出ASCII码（可以用于硬件如DSP、FPGA等）、.mat文件等等，这里就不再叙述。
利用MATLAB画出原始信号为

利用该滤波器对上述信号进行滤波可以得到如下

从图中可以看出，经过滤波器滤除之后的信号基本能够很好地恢复出原始的信号，从信号的频谱也能够看出，130Hz和1000Hz的信号被滤除。图中的频谱结果有一些幅度有一些衰减，这是正常的，同时频谱没有原始信号的频谱光滑，这是滤波器的阶数造成的，可以通过增加滤波器的阶数进行改善。
PS：由于滤波器存在一定的建立时间，即非理想滤波器，滤波之后的信号在刚开始的时间和后面的是不一致的，这里显示的是0.2s之后的信号。
代码如下：

clear;
close all;
clc;
f1=60;f2=130;f3=1e3;          %原始信号频率
fs=1e4;                       %采样频率
Ns=1e4;                       %采样点数
t=0:1/fs:Ns/fs;
s1=sin(2*pi*f1*t);
s2=sin(2*pi*f2*t);
s3=sin(2*pi*f3*t);
figure(1);
subplot(211)
plot(t,s1,'r',t,s2,'k',t,s3,'g');
xlim([0 500/fs])
title('正弦信号');xlabel('t/s');
y=s1+s2+s3;
set(gcf,'color','w');
subplot(212)
plot(t,y);
title('原始信号');xlabel('t/s');
set(gcf,'color','w');
N=4096;
f=(0:N-1)*fs/N;
yy=abs(fft(y,N)*2/N);
figure(2);
subplot(311)
plot(f(1:N/2),yy(1:N/2));
title('原始信号的FFT');xlabel('f/Hz');
xlim([0 1500]);
h=myfilter;
ylp=filter(h,y);
subplot(312)
plot(t,ylp);
xlim([0.2 0.8]);
title('低通滤波之后的信号')
xlabel('t/s');
yyy=abs(fft(ylp,N)*2/N);
subplot(313)
plot(f(1:N/2),yyy(1:N/2));
xlabel('f/Hz');
title('低通滤波之后的信号频谱');
xlim([0 1500]);