1. 从一次滑动平均优化说起认识filter函数的三种语法先讲个我自己的经历。早些年做传感器数据预处理任务很简单对一串带毛刺的电压信号做滑动平均。当时我第一反应是写个循环窗口长度设成10从第1个点滑到第N个点每个点取附近10个值的平均。循环写完后勉强能跑但数据一长就卡而且边界处理还得自己补逻辑。后来看到一个老工程师的代码里面就一行y filter(ones(1,10)/10, 1, x);我当时愣了一下随后查了文档才意识到MATLAB里专门干这件事的filter函数远比我想象的通用和强大。它不只是滑动平均而是一个能表达几乎所有线性时不变离散系统的通用工具。filter函数的作用用一句话说就是按照差分方程逐点计算输入信号通过一个离散系统后的输出。它本质上是求解线性常系数差分方程的工具也就是数字滤波器最底层的执行引擎。1.1 最常用的调用形态先看最常见的写法y filter(b, a, x)x输入信号通常是一个向量或矩阵b前馈系数分子多项式系数对应输入信号各延迟项的加权a反馈系数分母多项式系数对应输出信号历史值的加权y输出信号和x等长。它的数学本质是下面这个差分方程a(1)*y(n) b(1)*x(n) b(2)*x(n-1) ... b(nb1)*x(n-nb) - a(2)*y(n-1) - ... - a(na1)*y(n-na)如果看着有点晕换个角度理解filter做的事情就是对每一个新的输入点x(n)它用当前输入点和之前若干个输入点、之前若干个输出点的加权组合计算出新的输出点y(n)。这个加权组合的规律完全由b和a决定。拿滑动平均举例。5点滑动平均的公式是y(n) (x(n) x(n-1) x(n-2) x(n-3) x(n-4)) / 5写成filter就是b ones(1,5)/5; a 1; y filter(b, a, x);因为输出只依赖输入没有反馈项所以a就是1。这类滤波器叫有限脉冲响应滤波器FIR。一旦a里面除了1还有其他非零元素输出就会依赖之前的输出值这就成了无限脉冲响应滤波器IIR。1.2 带初始条件的调用形态很多人在初学阶段只用到上面的三参数形态但实际项目里尤其是实时处理场景三参数形态是不够的。原因在于filter默认把系统初始状态当成0也就是y(-1)0、y(-2)0……但真实系统在处理一段连续数据流时上一段末尾的状态会直接影响下一段开头的输出。这时候就要用四参数版本[y, zf] filter(b, a, x, zi)zi是初始状态zf是最终状态。这个机制我在后面坑的部分会详细展开这里先记住一个结论分段处理长信号、或者做实时流式滤波时必须用zi和zf衔接状态否则段与段的交界处会出现明显的跳变。还有一个用法是用filter画系统的脉冲响应和阶跃响应。本质上给filter输入一个单位脉冲序列[1, zeros(1,N)]输出就是系统的单位脉冲响应h(n)输入单位阶跃序列ones(1,N)输出就是阶跃响应。这个技巧在做系统分析时非常实用。N 50; imp [1, zeros(1, N-1)]; h filter(b, a, imp); % 脉冲响应 step ones(1, N); s filter(b, a, step); % 阶跃响应2. 先搞懂b和a再说调参一阶低通滤波器手把手拆解很多人用filter只会套滑动平均一旦遇到需要调滤波效果的场景就懵了。其实关键在于理解b和a的物理含义。我建议用一个最简单的一阶低通滤波器作为切入口把它彻底吃透其他滤波器自然就通了。2.1 一阶IIR滤波器的差分方程一个经典的一阶低通滤波器长这样y(n) alpha * y(n-1) (1 - alpha) * x(n)其中alpha在0到1之间。alpha越接近1滤波越平滑但响应越迟钝alpha越接近0滤波越弱但跟随越快。用filter实现这个公式时alpha 0.8; b 1 - alpha; % 也就是 0.2 a [1, -alpha]; % 也就是 [1, -0.8] y filter(b, a, x);为什么a是[1, -0.8]而不是[0.8]回到差分方程y(n) alpha * y(n-1) (1-alpha) * x(n)把y(n-1)项移到左边y(n) - alpha * y(n-1) (1-alpha) * x(n)对比filter的通用方程a(1)对应y(n)的系数1a(2)对应y(n-1)的系数就是-alpha。所以a [1, -alpha]注意第二项是负号。这个符号问题我见过太多人搞反写成了[1, alpha]结果滤波器直接发散。2.2 从Z变换视角理解系数如果想更深入一点可以看一眼Z变换。一阶低通滤波器的传递函数是H(z) (1 - alpha) / (1 - alpha * z^(-1))它在z alpha处有一个极点在z 0处有一个零点。当alpha从0向1变化时极点从原点向单位圆上的z1移动。极点在单位圆内系统稳定这对应alpha 1如果alpha 1极点跑到单位圆外滤波器就不稳定了输出会越来越大。用MATLAB可以直接可视化这个关系alpha 0.9; b 1 - alpha; a [1, -alpha]; zplane(b, a);zplane会画零极点图能非常直观地看到极点位置。这个图对理解滤波器稳定性特别有帮助。2.3 系数和截止频率的换算关系实际应用中我们更关心这个滤波器3dB截止频率是多少。一阶低通的alpha和截止频率fc归一化频率单位是cycles/sample之间有近似关系alpha exp(-2 * pi * fc)反过来fc -log(alpha) / (2 * pi)举个例子如果采样率是1000Hz想要截止频率在10Hz归一化截止频率就是10/1000 0.01那么fc_norm 0.01; alpha exp(-2*pi*fc_norm); % 约 0.9391 b 1 - alpha; a [1, -alpha];这样设计出来的滤波器大体上在10Hz附近开始衰减。需要注意的是这个公式只适用于一阶滤波器高阶滤波器的换算会复杂得多。遇到高阶需求直接用designfilt或fdatool设计系数再把生成的b、a丢给filter就行。3. filter的相关函数如何选择filter、filter2、filtfilt、conv的区别用filter一段时间后你一定会碰到filter2、filtfilt、conv、conv2这些长得相近的名字。它们到底有什么区别什么场景该用哪个我梳理了一张对比表然后逐个说明。函数适用维度核心用途输出长度是否因果filter1维每列独立按差分方程滤波IIR/FIR通吃与输入等长因果filter22维二维卷积滤波图像处理常用与输入等大因果默认filtfilt1维零相位滤波离线信号处理后处理与输入等长非因果conv1维计算两个序列的卷积length(x)length(h)-1非因果视h而定3.1 filter2是给图像准备的filter2和filter虽然只差一个数字但用途完全不同。filter2是二维滤波器专门处理矩阵信号最常见的应用就是图像滤波。比如给图像做一个3x3的均值模糊img imread(lena.png); img_gray rgb2gray(img); h ones(3,3)/9; img_blur filter2(h, img_gray, same);filter2的第三个参数控制输出尺寸same表示输出和输入一样大valid表示只输出卷积后完全重合的部分full表示完整卷积结果。图像处理里绝大多数场景用same就够了。需要注意的是filter2的输出是double类型如果后续要显示或转回uint8记得做类型转换和范围截取img_blur uint8(img_blur);3.2 filtfilt解决相位失真问题普通filter是因果滤波器每个时刻的输出只依赖当前及之前的输入和输出这必然带来相位延迟。信号中的特征点会被推后若干采样点。如果要精确分析波形的时间位置这个延时就很不友好。filtfilt的做法是先把信号正着过一遍滤波器再把结果倒过来过一遍同样的滤波器。两次滤波的相位延迟相互抵消输出信号就和原始信号的相位对齐了。y filtfilt(b, a, x);代价是filtfilt是非因果的它使用了未来的信息。所以它只能用于离线处理不能用于实时流式场景。还有一个使用前提信号长度要大于滤波器阶数的3倍否则两端的边界效应会非常明显。我之前处理心电信号时对比过用filter滤波后R波峰值位置偏移了约6个采样点换成filtfilt后偏差基本归零。如果你的场景是事后分析波形特征、测量时间间隔filtfilt是更好的选择如果场景是实时显示、在线控制那就只能用filter。3.3 conv和filter的微妙差异filter和conv都能实现FIR滤波但行为差别很大。conv(x, h)输出的是完整卷积长度是length(x)length(h)-1而filter(b, a, x)输出长度严格等于length(x)。如果你要做的只是FIR滤波且希望输出和输入等长大多数信号处理场景都是如此用filter。如果你需要计算两个信号的完整卷积结果或者想分析卷积的边界效应用conv。另一个区别是起始点的处理filter隐含了输入在n0之前都是0的假设输出从y(0)开始conv则从两个序列的第一个非零元素开始对齐计算。对同一个FIR滤波器系数做filter和conv前几个点的结果可能不一样这取决于你如何截取conv的结果。如果非要用conv模拟filter的效果需要自己把h补零对齐并截取中间段。我实际测试过用filter写起来更干净推荐直接使用filter。4. 项目实战filter函数处理三类常见信号任务第1到第3节把原理讲透了这一节直接上项目里能用的代码和思路。我挑三个最常见的场景实时去噪、差分方程求解、自定义FIR滤波器实现。4.1 场景一传感器数据的实时滑动平均去噪采集温度传感器数据时噪声通常是高频小幅波动滑动平均是最简单有效的去噪手段。但要注意如果你把整段数据一次性丢给filter那是离线处理真正的实时系统里数据是一个点一个点到达的必须用状态保持的方式。做法是先初始化状态然后每个新点到达时更新一次% 初始化 N 10; % 窗口长度 b ones(1, N)/N; a 1; zi zeros(size(b)-1); % FIR滤波器初始状态维度是 length(b)-1 % 每次来一个新数据点 for k 1:length(x_stream) [y_current, zi] filter(b, a, x_stream(k), zi); y_stream(k) y_current; end这里的关键是zi的类型。FIR滤波器的初始状态向量长度是length(b)-1IIR滤波器则是max(length(a), length(b))-1。搞错维度会直接报错正确写法可以这样统一生成zi zeros(max(length(a), length(b))-1, 1);如果不想手动管理状态MATLAB也提供了dsp.MovingAverage这类系统对象内部会维护状态。但从学习和理解角度我建议先用filter手动写一遍这样你对状态这个概念会有更扎实的认识。4.2 场景二用filter求解差分方程、模拟系统响应差分方程在数字信号处理、自动控制、经济学模型里到处都是。只要方程是线性、常系数的就能用filter求解。举个自动控制里的例子。设某离散系统的差分方程为y(n) - 0.6*y(n-1) 0.08*y(n-2) x(n) 0.5*x(n-1)要求系统在单位阶跃输入下的响应。直接编码b [1, 0.5]; a [1, -0.6, 0.08]; x ones(1, 100); % 单位阶跃100个点 y filter(b, a, x); plot(0:99, y); grid on;这个响应曲线就能告诉我们系统的稳态值、超调量、收敛速度等信息。稳态值可以直接算DC增益等于sum(b)/sum(a)也就是(10.5)/(1-0.60.08) 2.0833。用filter算出来的y(end)应该非常接近这个值。这个验证方法建议收藏任何时候不确定系数写没写对先算DC增益对比一下。4.3 场景三用filter实现自定义FIR滤波器如果你想设计一个带通滤波器比如保留500Hz到1000Hz之间的信号可以用designfilt先设计系数再交给filter执行。这是最稳妥的工程路径fs 8000; % 采样率8kHz % 设计带通滤波器阶数设为100 d designfilt(bandpassfir, FilterOrder, 100, ... CutoffFrequency1, 500, CutoffFrequency2, 1000, ... SampleRate, fs); b d.Coefficients; a 1; y filter(b, a, x);designfilt会返回一个滤波器对象其Coefficients属性就是FIR系数。对于FIR滤波器a永远是1所以filter(b, 1, x)就是完整的卷积滤波。我还想提醒一点用designfilt设计出的系数默认已经做了归一化DC增益是1。但如果你自己手动拼系数比如用窗函数法手工设计一定要检查增益。一个最简单的检查办法dc_gain sum(b) / sum(a); % 理想情况下等于1如果dc_gain偏离1很多信号经过滤波后整体就会被放大或缩小这在很多应用里是不允许的。5. 使用filter函数必须避开的四个坑最后这部分把我这些年用filter踩过的坑集中整理一下。这些坑在文档里不会写得特别起眼但在实际项目里都让我付出过代价。5.1 坑一a(1)必须为1否则行为容易混淆filter的文档规定a(1)不能为0但没说必须为1。实际上MATLAB会自动把整个差分方程除以a(1)来做归一化。这意味着filter([2, 4], [2, 6], x)和filter([1, 2], [1, 3], x)的结果是完全一样的。如果你在别的代码里看到a的第一个元素不是1不要以为它错了它只是没做归一化而已。但为了避免混淆我强烈建议所有代码一律写成a(1)1的形式。尤其是在不同代码块之间复制粘贴时两个写法看起来不一样容易让人误以为系数不同。5.2 坑二状态维度搞错导致报错或静默错误前面提过zi的维度问题。这里再补充一个更隐蔽的坑如果length(a)大于length(b)初始状态长度应该是length(a)-1反过来则是length(b)-1。统一公式是max(length(a), length(b)) - 1。我自己就写过这样的错误代码zi zeros(1, length(b)-1); % 假设 b 更长但某个滤波器实际是length(a)更长于是MATLAB报了维度不匹配。后来我改成统一计算zi zeros(max(length(a), length(b)) - 1, 1);就再也没在这个问题上翻过车。5.3 坑三边界瞬态被误判为真实信号filter默认初始状态为0这引入了一个从0开始的阶跃激励。如果滤波器阶数较高或者a系数的极点靠近单位圆输出开头一小段会经历明显的瞬态过渡幅值和相位都和稳态阶段不一样。这个瞬态不是输入信号的特征而是滤波器自身的启动过程。我在处理一段脉冲响应数据时差点把起始5个点的大幅波动当成了系统特性直到对比了filter和filtfilt的结果才意识到问题。解决方案有两个一是丢弃前若干点建议至少丢弃滤波器阶数个点二是用filtfilt做离线处理它会把边界效应集中到两端但仍然建议丢弃两端部分数据。5.4 坑四数据类型不注意精度悄悄丢失当x是single类型时filter返回的y也是single类型。这在信号较长、系数动态范围大的时候可能造成精度问题。另一个常见情况是x是整数类型比如uint8图像数据此时filter会强制把输出转成double但输入先被转成double再计算如果你的系数设计是基于double的结果没问题怕的是有人从某个CSV读入整数数据直接滤波后转回整数中间的截断误差被误认为是滤波效果。我的建议是进入filter之前统一把数据转成double计算完后再按需转回原类型。这种显式管理看起来多写两行代码但能省掉很多摸不着头脑的调试时间。5.5 补充一个实用技巧如何快速验证滤波器的正确性无论何时写完一个滤波器我建议立刻做两个验证第一构造一个单位脉冲输入算脉冲响应确认系统稳定响应趋向0第二构造一个直流信号全1算DC增益确认和解析值一致。两个都通过再去处理真实数据。这套流程每次大概花30秒但能拦住绝大部分低级错误。% 快速验证模板 imp [1, zeros(1, 99)]; step ones(1, 100); h filter(b, a, imp); s filter(b, a, step); assert(abs(sum(b)/sum(a) - s(end)) 1e-6, DC gain mismatch); if max(abs(h)) 1e6 warning(System seems unstable); endfilter这个函数本身很简单但它的强大在于能承载所有线性时不变系统的运算逻辑。从数字滤波器设计、系统仿真到信号预处理它都是底层引擎。搞懂它等于打通了MATLAB信号处理的一条主干道。
MATLAB filter函数完全指南：从滑动平均到IIR滤波的工程实践
1. 从一次滑动平均优化说起认识filter函数的三种语法先讲个我自己的经历。早些年做传感器数据预处理任务很简单对一串带毛刺的电压信号做滑动平均。当时我第一反应是写个循环窗口长度设成10从第1个点滑到第N个点每个点取附近10个值的平均。循环写完后勉强能跑但数据一长就卡而且边界处理还得自己补逻辑。后来看到一个老工程师的代码里面就一行y filter(ones(1,10)/10, 1, x);我当时愣了一下随后查了文档才意识到MATLAB里专门干这件事的filter函数远比我想象的通用和强大。它不只是滑动平均而是一个能表达几乎所有线性时不变离散系统的通用工具。filter函数的作用用一句话说就是按照差分方程逐点计算输入信号通过一个离散系统后的输出。它本质上是求解线性常系数差分方程的工具也就是数字滤波器最底层的执行引擎。1.1 最常用的调用形态先看最常见的写法y filter(b, a, x)x输入信号通常是一个向量或矩阵b前馈系数分子多项式系数对应输入信号各延迟项的加权a反馈系数分母多项式系数对应输出信号历史值的加权y输出信号和x等长。它的数学本质是下面这个差分方程a(1)*y(n) b(1)*x(n) b(2)*x(n-1) ... b(nb1)*x(n-nb) - a(2)*y(n-1) - ... - a(na1)*y(n-na)如果看着有点晕换个角度理解filter做的事情就是对每一个新的输入点x(n)它用当前输入点和之前若干个输入点、之前若干个输出点的加权组合计算出新的输出点y(n)。这个加权组合的规律完全由b和a决定。拿滑动平均举例。5点滑动平均的公式是y(n) (x(n) x(n-1) x(n-2) x(n-3) x(n-4)) / 5写成filter就是b ones(1,5)/5; a 1; y filter(b, a, x);因为输出只依赖输入没有反馈项所以a就是1。这类滤波器叫有限脉冲响应滤波器FIR。一旦a里面除了1还有其他非零元素输出就会依赖之前的输出值这就成了无限脉冲响应滤波器IIR。1.2 带初始条件的调用形态很多人在初学阶段只用到上面的三参数形态但实际项目里尤其是实时处理场景三参数形态是不够的。原因在于filter默认把系统初始状态当成0也就是y(-1)0、y(-2)0……但真实系统在处理一段连续数据流时上一段末尾的状态会直接影响下一段开头的输出。这时候就要用四参数版本[y, zf] filter(b, a, x, zi)zi是初始状态zf是最终状态。这个机制我在后面坑的部分会详细展开这里先记住一个结论分段处理长信号、或者做实时流式滤波时必须用zi和zf衔接状态否则段与段的交界处会出现明显的跳变。还有一个用法是用filter画系统的脉冲响应和阶跃响应。本质上给filter输入一个单位脉冲序列[1, zeros(1,N)]输出就是系统的单位脉冲响应h(n)输入单位阶跃序列ones(1,N)输出就是阶跃响应。这个技巧在做系统分析时非常实用。N 50; imp [1, zeros(1, N-1)]; h filter(b, a, imp); % 脉冲响应 step ones(1, N); s filter(b, a, step); % 阶跃响应2. 先搞懂b和a再说调参一阶低通滤波器手把手拆解很多人用filter只会套滑动平均一旦遇到需要调滤波效果的场景就懵了。其实关键在于理解b和a的物理含义。我建议用一个最简单的一阶低通滤波器作为切入口把它彻底吃透其他滤波器自然就通了。2.1 一阶IIR滤波器的差分方程一个经典的一阶低通滤波器长这样y(n) alpha * y(n-1) (1 - alpha) * x(n)其中alpha在0到1之间。alpha越接近1滤波越平滑但响应越迟钝alpha越接近0滤波越弱但跟随越快。用filter实现这个公式时alpha 0.8; b 1 - alpha; % 也就是 0.2 a [1, -alpha]; % 也就是 [1, -0.8] y filter(b, a, x);为什么a是[1, -0.8]而不是[0.8]回到差分方程y(n) alpha * y(n-1) (1-alpha) * x(n)把y(n-1)项移到左边y(n) - alpha * y(n-1) (1-alpha) * x(n)对比filter的通用方程a(1)对应y(n)的系数1a(2)对应y(n-1)的系数就是-alpha。所以a [1, -alpha]注意第二项是负号。这个符号问题我见过太多人搞反写成了[1, alpha]结果滤波器直接发散。2.2 从Z变换视角理解系数如果想更深入一点可以看一眼Z变换。一阶低通滤波器的传递函数是H(z) (1 - alpha) / (1 - alpha * z^(-1))它在z alpha处有一个极点在z 0处有一个零点。当alpha从0向1变化时极点从原点向单位圆上的z1移动。极点在单位圆内系统稳定这对应alpha 1如果alpha 1极点跑到单位圆外滤波器就不稳定了输出会越来越大。用MATLAB可以直接可视化这个关系alpha 0.9; b 1 - alpha; a [1, -alpha]; zplane(b, a);zplane会画零极点图能非常直观地看到极点位置。这个图对理解滤波器稳定性特别有帮助。2.3 系数和截止频率的换算关系实际应用中我们更关心这个滤波器3dB截止频率是多少。一阶低通的alpha和截止频率fc归一化频率单位是cycles/sample之间有近似关系alpha exp(-2 * pi * fc)反过来fc -log(alpha) / (2 * pi)举个例子如果采样率是1000Hz想要截止频率在10Hz归一化截止频率就是10/1000 0.01那么fc_norm 0.01; alpha exp(-2*pi*fc_norm); % 约 0.9391 b 1 - alpha; a [1, -alpha];这样设计出来的滤波器大体上在10Hz附近开始衰减。需要注意的是这个公式只适用于一阶滤波器高阶滤波器的换算会复杂得多。遇到高阶需求直接用designfilt或fdatool设计系数再把生成的b、a丢给filter就行。3. filter的相关函数如何选择filter、filter2、filtfilt、conv的区别用filter一段时间后你一定会碰到filter2、filtfilt、conv、conv2这些长得相近的名字。它们到底有什么区别什么场景该用哪个我梳理了一张对比表然后逐个说明。函数适用维度核心用途输出长度是否因果filter1维每列独立按差分方程滤波IIR/FIR通吃与输入等长因果filter22维二维卷积滤波图像处理常用与输入等大因果默认filtfilt1维零相位滤波离线信号处理后处理与输入等长非因果conv1维计算两个序列的卷积length(x)length(h)-1非因果视h而定3.1 filter2是给图像准备的filter2和filter虽然只差一个数字但用途完全不同。filter2是二维滤波器专门处理矩阵信号最常见的应用就是图像滤波。比如给图像做一个3x3的均值模糊img imread(lena.png); img_gray rgb2gray(img); h ones(3,3)/9; img_blur filter2(h, img_gray, same);filter2的第三个参数控制输出尺寸same表示输出和输入一样大valid表示只输出卷积后完全重合的部分full表示完整卷积结果。图像处理里绝大多数场景用same就够了。需要注意的是filter2的输出是double类型如果后续要显示或转回uint8记得做类型转换和范围截取img_blur uint8(img_blur);3.2 filtfilt解决相位失真问题普通filter是因果滤波器每个时刻的输出只依赖当前及之前的输入和输出这必然带来相位延迟。信号中的特征点会被推后若干采样点。如果要精确分析波形的时间位置这个延时就很不友好。filtfilt的做法是先把信号正着过一遍滤波器再把结果倒过来过一遍同样的滤波器。两次滤波的相位延迟相互抵消输出信号就和原始信号的相位对齐了。y filtfilt(b, a, x);代价是filtfilt是非因果的它使用了未来的信息。所以它只能用于离线处理不能用于实时流式场景。还有一个使用前提信号长度要大于滤波器阶数的3倍否则两端的边界效应会非常明显。我之前处理心电信号时对比过用filter滤波后R波峰值位置偏移了约6个采样点换成filtfilt后偏差基本归零。如果你的场景是事后分析波形特征、测量时间间隔filtfilt是更好的选择如果场景是实时显示、在线控制那就只能用filter。3.3 conv和filter的微妙差异filter和conv都能实现FIR滤波但行为差别很大。conv(x, h)输出的是完整卷积长度是length(x)length(h)-1而filter(b, a, x)输出长度严格等于length(x)。如果你要做的只是FIR滤波且希望输出和输入等长大多数信号处理场景都是如此用filter。如果你需要计算两个信号的完整卷积结果或者想分析卷积的边界效应用conv。另一个区别是起始点的处理filter隐含了输入在n0之前都是0的假设输出从y(0)开始conv则从两个序列的第一个非零元素开始对齐计算。对同一个FIR滤波器系数做filter和conv前几个点的结果可能不一样这取决于你如何截取conv的结果。如果非要用conv模拟filter的效果需要自己把h补零对齐并截取中间段。我实际测试过用filter写起来更干净推荐直接使用filter。4. 项目实战filter函数处理三类常见信号任务第1到第3节把原理讲透了这一节直接上项目里能用的代码和思路。我挑三个最常见的场景实时去噪、差分方程求解、自定义FIR滤波器实现。4.1 场景一传感器数据的实时滑动平均去噪采集温度传感器数据时噪声通常是高频小幅波动滑动平均是最简单有效的去噪手段。但要注意如果你把整段数据一次性丢给filter那是离线处理真正的实时系统里数据是一个点一个点到达的必须用状态保持的方式。做法是先初始化状态然后每个新点到达时更新一次% 初始化 N 10; % 窗口长度 b ones(1, N)/N; a 1; zi zeros(size(b)-1); % FIR滤波器初始状态维度是 length(b)-1 % 每次来一个新数据点 for k 1:length(x_stream) [y_current, zi] filter(b, a, x_stream(k), zi); y_stream(k) y_current; end这里的关键是zi的类型。FIR滤波器的初始状态向量长度是length(b)-1IIR滤波器则是max(length(a), length(b))-1。搞错维度会直接报错正确写法可以这样统一生成zi zeros(max(length(a), length(b))-1, 1);如果不想手动管理状态MATLAB也提供了dsp.MovingAverage这类系统对象内部会维护状态。但从学习和理解角度我建议先用filter手动写一遍这样你对状态这个概念会有更扎实的认识。4.2 场景二用filter求解差分方程、模拟系统响应差分方程在数字信号处理、自动控制、经济学模型里到处都是。只要方程是线性、常系数的就能用filter求解。举个自动控制里的例子。设某离散系统的差分方程为y(n) - 0.6*y(n-1) 0.08*y(n-2) x(n) 0.5*x(n-1)要求系统在单位阶跃输入下的响应。直接编码b [1, 0.5]; a [1, -0.6, 0.08]; x ones(1, 100); % 单位阶跃100个点 y filter(b, a, x); plot(0:99, y); grid on;这个响应曲线就能告诉我们系统的稳态值、超调量、收敛速度等信息。稳态值可以直接算DC增益等于sum(b)/sum(a)也就是(10.5)/(1-0.60.08) 2.0833。用filter算出来的y(end)应该非常接近这个值。这个验证方法建议收藏任何时候不确定系数写没写对先算DC增益对比一下。4.3 场景三用filter实现自定义FIR滤波器如果你想设计一个带通滤波器比如保留500Hz到1000Hz之间的信号可以用designfilt先设计系数再交给filter执行。这是最稳妥的工程路径fs 8000; % 采样率8kHz % 设计带通滤波器阶数设为100 d designfilt(bandpassfir, FilterOrder, 100, ... CutoffFrequency1, 500, CutoffFrequency2, 1000, ... SampleRate, fs); b d.Coefficients; a 1; y filter(b, a, x);designfilt会返回一个滤波器对象其Coefficients属性就是FIR系数。对于FIR滤波器a永远是1所以filter(b, 1, x)就是完整的卷积滤波。我还想提醒一点用designfilt设计出的系数默认已经做了归一化DC增益是1。但如果你自己手动拼系数比如用窗函数法手工设计一定要检查增益。一个最简单的检查办法dc_gain sum(b) / sum(a); % 理想情况下等于1如果dc_gain偏离1很多信号经过滤波后整体就会被放大或缩小这在很多应用里是不允许的。5. 使用filter函数必须避开的四个坑最后这部分把我这些年用filter踩过的坑集中整理一下。这些坑在文档里不会写得特别起眼但在实际项目里都让我付出过代价。5.1 坑一a(1)必须为1否则行为容易混淆filter的文档规定a(1)不能为0但没说必须为1。实际上MATLAB会自动把整个差分方程除以a(1)来做归一化。这意味着filter([2, 4], [2, 6], x)和filter([1, 2], [1, 3], x)的结果是完全一样的。如果你在别的代码里看到a的第一个元素不是1不要以为它错了它只是没做归一化而已。但为了避免混淆我强烈建议所有代码一律写成a(1)1的形式。尤其是在不同代码块之间复制粘贴时两个写法看起来不一样容易让人误以为系数不同。5.2 坑二状态维度搞错导致报错或静默错误前面提过zi的维度问题。这里再补充一个更隐蔽的坑如果length(a)大于length(b)初始状态长度应该是length(a)-1反过来则是length(b)-1。统一公式是max(length(a), length(b)) - 1。我自己就写过这样的错误代码zi zeros(1, length(b)-1); % 假设 b 更长但某个滤波器实际是length(a)更长于是MATLAB报了维度不匹配。后来我改成统一计算zi zeros(max(length(a), length(b)) - 1, 1);就再也没在这个问题上翻过车。5.3 坑三边界瞬态被误判为真实信号filter默认初始状态为0这引入了一个从0开始的阶跃激励。如果滤波器阶数较高或者a系数的极点靠近单位圆输出开头一小段会经历明显的瞬态过渡幅值和相位都和稳态阶段不一样。这个瞬态不是输入信号的特征而是滤波器自身的启动过程。我在处理一段脉冲响应数据时差点把起始5个点的大幅波动当成了系统特性直到对比了filter和filtfilt的结果才意识到问题。解决方案有两个一是丢弃前若干点建议至少丢弃滤波器阶数个点二是用filtfilt做离线处理它会把边界效应集中到两端但仍然建议丢弃两端部分数据。5.4 坑四数据类型不注意精度悄悄丢失当x是single类型时filter返回的y也是single类型。这在信号较长、系数动态范围大的时候可能造成精度问题。另一个常见情况是x是整数类型比如uint8图像数据此时filter会强制把输出转成double但输入先被转成double再计算如果你的系数设计是基于double的结果没问题怕的是有人从某个CSV读入整数数据直接滤波后转回整数中间的截断误差被误认为是滤波效果。我的建议是进入filter之前统一把数据转成double计算完后再按需转回原类型。这种显式管理看起来多写两行代码但能省掉很多摸不着头脑的调试时间。5.5 补充一个实用技巧如何快速验证滤波器的正确性无论何时写完一个滤波器我建议立刻做两个验证第一构造一个单位脉冲输入算脉冲响应确认系统稳定响应趋向0第二构造一个直流信号全1算DC增益确认和解析值一致。两个都通过再去处理真实数据。这套流程每次大概花30秒但能拦住绝大部分低级错误。% 快速验证模板 imp [1, zeros(1, 99)]; step ones(1, 100); h filter(b, a, imp); s filter(b, a, step); assert(abs(sum(b)/sum(a) - s(end)) 1e-6, DC gain mismatch); if max(abs(h)) 1e6 warning(System seems unstable); endfilter这个函数本身很简单但它的强大在于能承载所有线性时不变系统的运算逻辑。从数字滤波器设计、系统仿真到信号预处理它都是底层引擎。搞懂它等于打通了MATLAB信号处理的一条主干道。
