资讯动态

MATLAB 5 点滑动平均:filter 与 conv 的输出对齐,以及 2 样本延时的含义

发布时间:2026/10/14 11:59:21 来源:尧图企业网站定制
MATLAB 5 点滑动平均filter 与 conv 的输出对齐以及 2 样本延时的含义filter(ones(1,5)/5,1,x)和conv(x,ones(1,5)/5,same)都能得到一条平滑曲线长度也都与 x 相同但同一数组下标上的结果不一定相等。另一方面5 点对称 FIR 的群延时是 2 个样本也不代表输出就是把输入原样向后挪两格。这篇用冲激、阶跃和两个频率叠加的信号分别检查输出截取、启动暂态和频率相关增益。代码只用基础 MATLAB不调用freqz、grpdelay或filtfilt。1. 运行方式与实测结果把第 3 节完整代码保存为example03_average.m在 MATLAB 中切换到该文件所在目录后执行example03_average终端配置好 MATLAB 命令行路径后也可执行matlab -batch example03_average本次实际测试版本为25.2.0.2998904 (R2025b)。所有输入在函数内构造没有外部数据。正常运行会在后台导出output03/moving_average.png并保存数值指标到output03/average_metrics.csv。同一输出目录重跑会覆盖示例文件。需要保留不同运行时可以执行example03_average(output03_run2,false)。2026 年 10 月 14 日正向 batch 退出码为 010 项断言通过。关键结果如下。检查项实测结果filter 与 full 卷积前 N 项的最大绝对误差2.22044604925e-16N100 时 full 卷积长度104单位阶跃前五个输出0.2、0.4、0.6、0.8、1稳态两频信号与解析式的最大绝对误差4.60742555219e-15稳态输出与单位增益纯延时的最大差异0.8103722208310.08π rad/sample 对应的有符号增益0.9379559364690.6π rad/sample 对应的有符号增益-0.2472135955避开频响零点后计算的群延时范围1.99999999996 至 2.00000000001 样本很小的数值误差反映 double 浮点运算的舍入0.81037 的差异则说明“只做两样本纯延时”的模型与平均滤波器并不相同。2. 先说明时间坐标和边界五点因果平均定义为y[n] (x[n] x[n-1] x[n-2] x[n-3] x[n-4])/5。每个输出使用当前样本和前四个样本。本文令记录前的输入为零采用filter的默认零初始状态。这种初始化和输出同长规则见 filter R2025b 文档。长度 N 的 x 与长度 5 的 b 做 full 卷积共有 N4 项。filter(b,1,x)输出前 N 项full 卷积还保留输入记录结束后、历史数据逐步离开窗口的四项尾部。same取的是 full 卷积的中心 N 项。对本例的五点核相当于 full 结果的第 3 到第 N2 项而因果记录取第 1 到第 N 项。因此把两者都画在原来的n0:N-1上会产生两样本的对齐差异末尾处理也不相同。same的定义可查 conv R2025b 文档。如果把居中结果放在原输入时刻上理解它需要使用该时刻之后的两个样本用于在线处理时需要明确等待这些样本的时间。3. 完整可运行代码以下为实际测试文件。x 和 b 均为行向量避免比较时混入行列方向问题。functionexample03_average(outDir,negativeMode)% Causal 5-point moving average, zero initial state. Base MATLAB only.ifnargin1,outDirfullfile(fileparts(mfilename(fullpath)),output03);endifnargin2,negativeModefalse;endif~isfolder(outDir),mkdir(outDir);endfprintf(MATLAB_VERSION%s\nNEGATIVE_MODE%d\n,version,negativeMode);tol1e-10;bones(1,5)/5;N100;n0:N-1;w10.08*pi;w20.6*pi;xcos(w1*n)0.6*cos(w2*n);yfilter(b,1,x);fullYconv(x,b,full);pureDelay[00x(1:end-2)];ifnegativeMode,ypureDelay;end% Intentionally incorrect only a delay.convErrormax(abs(y-fullY(1:N)));ck(convErrortol,Convolution,Causal filter equals first N terms of full convolution);ck(numel(fullY)N4,FullLength,Full convolution includes the four-sample tail);sameYconv(x,b,same);ck(max(abs(y-sameY))0.1,NotSame,Centered same convolution is not the causal record);impulsefilter(b,1,[1zeros(1,11)]);ck(max(abs(impulse-[bzeros(1,7)]))tol,Impulse,Impulse response has five taps of 0.2);stepfilter(b,1,ones(1,12));ck(max(abs(step(1:5)-(1:5)/5))tol,Transient,Zero-state step front is a ramp);ck(max(abs(step(5:end)-1))tolabs(sum(b)-1)tol,DC,Steady-state DC gain is one);% H(w)exp(-j*2*w)*A(w), A real, possibly negative.gain(w)sin(5*w/2)./(5*sin(w/2));predictedgain(w1)*cos(w1*(n-2))0.6*gain(w2)*cos(w2*(n-2));steadyErrormax(abs(y(5:end)-predicted(5:end)));ck(steadyErrortol,SteadyTone,Steady sinusoids have frequency-dependent signed gain);delayErrormax(abs(y(5:end)-pureDelay(5:end)));ck(delayError0.1,NotPureDelay,Mixed frequencies are not a unit-gain pure shift);% Compute phase derivative analytically, without freqz or grpdelay.omegalinspace(0.001,pi,2000);Eexp(-1i*omega*(0:4));HE*b;Hprime-1i*E*((0:4).*b);awayabs(H)1e-3;gd-imag(Hprime(away)./H(away));ck(max(abs(gd-2))1e-9,GroupDelay,Group delay is two away from response zeros);zeroHexp(-1i*[2*pi/5;4*pi/5]*(0:4))*b;ck(max(abs(zeroH))tol,Zeros,Response zeros: group delay is undefined there);fprintf(CONV_MAX_ERROR%.12g STEADY_MAX_ERROR%.12g PURE_DELAY_MAX_ERROR%.12g\n,convError,steadyError,delayError);fprintf(GAIN_LOW%.12g GAIN_HIGH%.12g GROUP_DELAY_RANGE[%.12g %.12g]\n,gain(w1),gain(w2),min(gd),max(gd));fprintf(STEP_FRONT%s\nIMPULSE%s\n,mat2str(step(1:5),8),mat2str(impulse,8));fprintf(LIMIT: group delay only where H is nonzero; startup and frequency-dependent gains prevent generic pure shifts.\n);figfigure(Visible,off,Color,w,Position,[1001001100760]);cleanuponCleanup(()close(fig));tiledlayout(2,2,TileSpacing,compact);nexttile;stem(0:11,impulse,filled);grid on;ylim([00.25]);title(Impulse: five coefficients of 0.2);xlabel(Sample n);ylabel(h[n]);nexttile;stairs(0:11,step,LineWidth,1.5);grid on;ylim([01.1]);title(Zero-state step: startup ramp then unity);xlabel(Sample n);ylabel(Output);nexttile;plot(n,y,n,pureDelay,--,LineWidth,1.1);xlim([045]);grid on;legend(5-point average,Pure delay of 2,Location,best);title(Same nominal delay does not mean same waveform);xlabel(Sample n);ylabel(Amplitude);nexttile;plot(omega/pi,abs(H),LineWidth,1.3);grid on;xlabel(Normalized frequency (x pi rad/sample));ylabel(|H|);title(Frequency-dependent gain; zeros at 0.4 and 0.8);exportgraphics(fig,fullfile(outDir,moving_average.png),Resolution,150);metricstable(convError,steadyError,delayError,min(gd),max(gd),...VariableNames,{conv_error,steady_error,pure_delay_error,gd_min,gd_max});writetable(metrics,fullfile(outDir,average_metrics.csv));fprintf(RESULTPASS example03_average\n);endfunctionck(ok,id,msg)assert(ok,[lessonMA:id],msg);fprintf(ASSERT_PASS%s\n,id);end4. 从冲激和阶跃读懂初始状态左上图中单位冲激只有第一个输入为 1后面全部为 0。这个 1 在五个连续窗口中各贡献 1/5因此冲激响应前五项都是 0.2之后为 0。这也直接对应代码里的五个滤波系数。右上图输入是从记录开始就为 1 的单位阶跃。由于记录前四个输入按零处理第一项是(10000)/50.2之后依次为 0.4、0.6、0.8到第五个输出才达到 1。图中的横轴 n 从 0 开始所以第五个输出位于 n4。它与 MATLAB 数组里的step(5)是同一个点别在读图时多算一拍。系数和为 1说明稳态直流增益为 1启动阶段没有填满窗口仍会有这段暂态。如果希望从第一项就符合另一种历史输入需要明确初始条件不能只把曲线开头的偏低判成算法错误。默认movmean(x,5)也不等价于这里的零初态因果平均它采用居中窗口端点默认缩小窗口。movmean 的端点说明可用于核对这类差别。5. 为什么群延时是 2波形却没有原样平移五点平均的频率响应可写为H(ω) exp(-j*2*ω) * sin(5*ω/2)/(5*sin(ω/2))。其中 ω 的单位是 rad/sampleω0 时按极限得到 H(0)1。前面的指数项对应两样本延时后面的实因子随频率变化并且可能为负。在 H 不为零的连续频段上这个实因子的符号不会突然改变相位对 ω 的斜率为 -2因此群延时为 2。跨过零点会发生相位跳变零点处输出分量被消除群延时没有可用定义。群延时与相位延时的概念可参阅 Group Delay and Phase Delay。本例不用数值差分相位而是直接对有限 FIR 和求导计算 H、H′再由-imag(Hprime./H)得到群延时。为了避开零点附近的除法放大检查范围限定为abs(H)1e-3。实测范围非常接近 2同时另有断言检查 0.4π、0.8π 两个零点。5.1 两个频率会受到不同衰减输入由下面两部分相加而成低频余弦ω₁0.08π输入幅度为 1高频余弦ω₂0.6π输入幅度为 0.6进入稳态后低频项的有符号增益约为 0.937956高频项约为 -0.247214。负号表示相对于共同的两样本延时形式还存在符号反转它不代表幅频响应为负。右下图画的是abs(H)因此始终非负。两个分量经历不同的衰减和符号变化叠加出来的形状自然会改变。左下图蓝线是平均滤波器输出虚线是把原始输入仅向后移两格两条曲线明显不同。代码从第五项开始比较稳态解析式最大误差约为 4.61e-15与纯延时的最大差异约为 0.81037。因此群延时描述的是相位随频率变化的性质。只有在相关频带的幅度响应也足够接近所需增益、并妥善处理暂态时才能把某些输入近似理解为延后而形状基本保持。5.2 换算成秒时要用采样率两样本的算法群延时对应2/Fs秒。假设采样率为 1 kHz这个换算值是 2 ms这是单位换算举例本文代码没有指定 1 kHz 的硬件采样率。如果实现到 FPGA额外寄存器流水线通常还会引入系统时钟周期延迟。采样间隔与系统时钟周期可能不同不能把“2 个样本”直接说成“2 个 FPGA 时钟周期”。6. 让“只是延时”的错误模型失败执行example03_average(output03_negative,true)负对照把 y 故意替换为两样本纯延时再与 full 卷积的前 N 项比较。本次实际在lessonMA:Convolution失败batch 退出码为 1。这个负对照验证的是错误的数学模型能否被检查抓住它没有运行硬件也没有证明任何固定点位宽、溢出策略或 RTL 实现正确。7. 使用范围和参考资料本例只验证了零初态、固定五点系数、double 浮点输入。平均会削弱快速变化也会让边沿展开。不能只凭平滑后的曲线更顺眼就认为有效信号全部保留了。分块处理时需要传递滤波器状态否则每块开头都会重新产生零初态暂态。官方filter文档给出了zi/zf用法本文代码没有运行分块状态实验也没有把它列入上述 10 项断言。选题参考 Convolution in Digital Signal Processing 固定版本 README中的卷积、LTI 与平均滤波概念。完整课程使用 Signal Processing Toolbox部分案例另需图像处理和深度学习工具箱。该版本同时包含 BSD 3-Clause 文本的 LICENSE.md以及带 MathWorks 产品及服务使用限制的 license.txt。本文仅作概念参考实验代码和图片独立生成不复制仓库素材。需要进一步区分算法延时和硬件流水线延时可以阅读 流水线设计 Pipeline提升时序性能的核心技巧。最后留一个具体问题如果把同一段输入分成两块第二块仍直接调用没有初始状态参数的filter它的前四个输出为什么可能与整段处理不同这与本篇阶跃前沿的四个过渡点有什么联系

读完文章,也想定制专属网站?

尧图设计师 24 小时内与您沟通定制方案

免费获取报价 →
↑