资讯动态

MATLAB 分块 FIR 滤波如何保持连续:用 zi/zf 消除重复启动瞬态

发布时间:2026/10/18 3:39:24 来源:尧图企业网站定制
实测说明2026-10-17本文代码已在 Windows 11构建 26200、MATLAB 25.2.0.2998904R2025b、PCWIN64 中实际执行主实验及独立边界补测均通过。下列表格和配图来自真实 MATLAB 输出requiredFilesAndProducts 实际分析仅返回基础 MATLAB不需要 Signal Processing Toolbox。另保留 Python 3.12.14、NumPy 2.3.5、SciPy 1.17.0 的先期参考检查作为独立补充证据。最短运行将后文完整代码保存为block_fir_state_demo.m在文件所在目录执行results block_fir_state_demo;。真实关键结果保留状态的分块输出与整段输出最大差为 2.22×10⁻¹⁶每块清零时最大误差为 1.149161978。从 FPGA 经 DMA、串口或网络收到一块数据就调用一次filter(b,1,xBlock)很容易得到一种迷惑性的结果块中间看起来正常块头总有一段异常。改变块长以后异常出现的位置也跟着变化。先检查每次调用是否保留了滤波器状态。对同一条连续采样流软件收包边界通常没有物理含义。如果每块都从零状态开始算法会把每个包都当成一次新的开机采样。1. 一条记录切成几块应该得到同一组输出本例固定采样率Fs2000 Hz记录长度N4097使用 33 抽头对称 FIR。系数由理想低通冲激响应乘 Hamming 形式的窗口后归一化生成全部按公式编写不调用fir1、hamming或sinc。这里的 180 Hz 是构造理想低通响应时的截止参数。有限长度加窗后实际通带纹波、过渡带和阻带衰减要另行测量本例不把它当作经过指标设计的产品滤波器。输入包含 41 Hz 主信号、647 Hz 分量、一次直流阶跃和一个孤立脉冲。分块长度循环取pattern[176432571289];这组长度故意包含单点块、比 FIR 记忆长度更短的块以及不能整除总记录长度的块。共得到 62 块。它比只测试“每块 1024 点”更容易暴露状态和最后一块的处理问题。我们比较三条路径整段输入一次filter分块输入把上一块zf传给下一块zi分块输入每次清零状态作为故意错误的对照前两条输出应在浮点误差范围内相同。第三条在块头反复丢失历史信息。2. zi 保存的是滤波器内部状态FIR 的因果关系为y [ n ] ∑ k 0 L − 1 b [ k ] x [ n − k ] y[n]\sum_{k0}^{L-1}b[k]x[n-k]y[n]k0∑L−1​b[k]x[n−k]当L33时当前输出可能依赖前 32 个输入样本。块头没有这些历史输入直接重新计算就等价于把它们设为零。MATLAB 提供的接口为[yBlock,zf]filter(b,1,xBlock,zi,1);zizf;对这里的 FIR状态长度为L-132。MathWorks 文档明确说明了初始状态、最终状态和沿指定维度过滤的用法。MathWorks filter 文档需要留意zi是filter所用实现结构的延迟状态不能直接用“最近 32 个原始采样值”替换。若在自己的 FPGA FIR 中保存输入移位寄存器那是另一种状态表示要么比较两边的输出要么推导结构之间的状态映射。独立补测还在 MATLAB 中手写转置直接 II 型 FIR 状态递推逐样本核对 filter 输出输出最大差为 4.44×10⁻¹⁶末状态最大差为 1.11×10⁻¹⁶。先期 Python 参考模型也做过同类核对两种环境的记录分别保存未把 Python 输出当作 MATLAB 实测。3. 多通道时单点块最容易触发维度问题约定数据始终为N×C每行一个时刻每列一个通道。两路信号对应两套独立状态zizeros(numel(b)-1,C);[yBlock,zi]filter(b,1,xBlock,zi,1);最后的1明确指定沿时间维处理。特别是xBlock只有一行时它的形状是1×C不能再依赖“第一个非单一维度”的默认规则。否则本应彼此独立的通道有可能被当成同一条序列的相邻样本。同样不能为了修复维度直接对多通道矩阵调用xBlock(:)。这种写法会改变样本与通道的组织关系。单通道行向量可以在输入接口处转成列向量多通道数据则必须明确原始布局。4. 16 点群延迟和 32 点历史记忆是两件事这组系数对称长度为 33因此其线性相位群延迟为(33-1)/216个采样点也就是 8 ms。该结论用于滤波器响应不为零、相位有定义的频段。为了把有限记录的完整线性卷积输出取完需要在真实流结束后继续送入L-132个零。16 和 32 分别对应群延迟与冲激响应尾部长度不能互相替代。代码先比较N个因果输出再用保存的状态处理 32 个零[yTail,zEnd]filter(b,1,zeros(L-1,1),zf,1);yFull[yStream;yTail];yFull应与conv(x,b)的完整结果一致。conv的完整输出长度为两个输入长度之和减一。MathWorks conv 文档连续采样过程中不要每收完一个包就冲刷 32 个零。那相当于人为往物理时间线上插入零输入。只有真正结束一条有限流、并且明确要输出尾部时才做这一步。5. 本次 MATLAB 实测结果检查项本次 MATLAB 结果变长块保留状态与整段滤波最大差值2.22×10⁻¹⁶保留状态与直接卷积前 N 点最大差值8.88×10⁻¹⁶每块清零与整段滤波最大差值1.149161978分块与整段滤波的最终状态差值0最后一次零冲刷后的完整卷积差值8.88×10⁻¹⁶多通道结果最大差值2.22×10⁻¹⁶六种不同分块方式中的最大差值2.22×10⁻¹⁶独立补测手写逐点递推与 filter 最大差值4.44×10⁻¹⁶冲激输入与 FIR 系数序列差值0MATLAB 独立补测验证了长度为 1、2、7、31、32、33、34 的记录每种长度分别输入零、常量和确定性信号共 21 组每组分别检查整段输出、末状态、完整卷积尾部和补尾后状态。还测试了每次只输入一个样本的双通道记录以及错误清零方法的指定反例。上述检查全部通过。本文主程序的 9 处断言均实际执行一次。与同批联合拟合实验合计两个主程序有 20 次成功断言执行独立补测有 103 个通过检查这些计数不是 103 个独立数据集也不包括额外的两次 MAT 保存回读核验。浮点运算顺序不同末位误差可以变化本例采用1e-12作为该组双精度数据的数值容差并不把它作为任意幅度、任意长度数据的通用标准。图中灰色竖线是分块起点。第一幅可以看到错误输出在部分边界突然回到零附近第二幅直接画误差第三幅验证整条有限流结束后的输出尾部。它们均为本次 MATLAB 实际导出的图。6. 完整 MATLAB 示例与运行方法保存下面代码为block_fir_state_demo.m在文件所在目录运行resultsblock_fir_state_demo;函数会创建matlab_results目录执行断言并输出 PNG、MAT 文件。本次实际日志已打印BLOCK MATLAB PASS; blocks62; taps33; group_delay_samples16主程序批次进程退出码为 0。复现时可先运行version、ver和which filter -all记录环境排查同名函数遮蔽使用新的输出目录保存每次结果。functionresultsblock_fir_state_demo(outDir)% Base-MATLAB FIR streaming demo with independent channel states.% Compares retained states, reset states, partitions, and finite-stream tails.ifnargin1outDirfullfile(fileparts(mfilename(fullpath)),matlab_results);endif~exist(outDir,dir),mkdir(outDir);endFs2000;N4097;L33;fc180;n(0:N-1);tn/Fs;bmake_lowpass(L,fc/Fs);x0.8*sin(2*pi*41*t0.2)0.2*cos(2*pi*647*t);x(n1550)x(n1550)0.35;x(2501)x(2501)1.1;pattern[176432571289];[yWhole,zWhole]filter(b,1,x,zeros(L-1,1),1);[yKeep,zKeep,starts]run_blocks(x,b,pattern,true);[yReset,~]run_blocks(x,b,pattern,false);yConvconv(x,b(:));keepErrormax(abs(yKeep-yWhole));convErrormax(abs(yKeep-yConv(1:N)));resetErrormax(abs(yReset-yWhole));stateErrormax(abs(zKeep-zWhole));assert(keepError1e-12,BLOCK:Continuity,Block continuity failed.);assert(convError1e-12,BLOCK:Convolution,Causal convolution differs.);assert(stateError1e-12,BLOCK:State,Final states differ.);assert(resetError0.1,BLOCK:NegativeControl,Reset error not exposed.);% Flush exactly L-1 zero inputs only when the finite stream really ends.[yTail,zAfterTail]filter(b,1,zeros(L-1,1),zKeep,1);yFull[yKeep;yTail];tailErrormax(abs(yFull-yConv));assert(tailError1e-12,BLOCK:Tail,Full convolution tail differs.);assert(max(abs(zAfterTail))1e-12,BLOCK:Flush,State did not clear.);% Every column is a separate channel. dim1 also covers 1-by-C blocks.x2[x,0.4*cos(2*pi*73*t)-0.15];y2Wholefilter(b,1,x2,zeros(L-1,2),1);[y2Keep,~]run_blocks(x2,b,pattern,true);channelErrormax(abs(y2Keep(:)-y2Whole(:)));assert(channelError1e-12,BLOCK:Channels,Channel states mixed.);% Same input, different partitions; include one sample per call.patterns{1,[235],64,1024,[313233],N};partitionError0;fork1:numel(patterns)yTestrun_blocks(x,b,patterns{k},true);partitionErrormax(partitionError,max(abs(yTest-yWhole)));endassert(partitionError1e-12,BLOCK:Partitions,Partition changed output.);% An impulse at sample zero must yield the tap sequence, followed by zeros.xImpulse[1;zeros(2*L-1,1)];yImpulserun_blocks(xImpulse,b,[1311],true);impulseExpected[b(:);zeros(L,1)];impulseErrormax(abs(yImpulse-impulseExpected));assert(impulseError1e-12,BLOCK:Impulse,Impulse check failed.);figfigure(Visible,off,Color,w,Position,[1001001050760]);subplot(3,1,1);plot(n,yWhole,k-,n,yReset,r-);grid on;xlim([0650]);ylabel(Output);legend(Continuous state,Reset every block);title(FIR block boundaries: same samples, different state policy);subplot(3,1,2);plot(n,yReset-yWhole,r-);hold on;plot(n,yKeep-yWhole,b-);grid on;xlim([0650]);ylabel(Error);legend(Reset - whole,Keep - whole);ylim([-1.41.4]);fork1:numel(starts)ifstarts(k)-1650,break;endplot([starts(k)-1starts(k)-1],[-1.41.4],:,...Color,[0.70.70.7],HandleVisibility,off);endsubplot(3,1,3);plot(0:numel(yFull)-1,yConv,k-,0:numel(yFull)-1,yFull,b--);xlim([N-45NL-2]);grid on;xlabel(Zero-based sample index);ylabel(Output);legend(Full convolution,Stream final zero flush);title(Only append the FIR tail at the true end of a finite record);print(fig,fullfile(outDir,block_fir_state_matlab.png),-dpng,-r150);close(fig);resultsstruct(Fs,Fs,N,N,taps,L,groupDelaySamples,(L-1)/2,...blocks,numel(starts),keepError,keepError,convError,convError,...resetError,resetError,stateError,stateError,tailError,tailError,...channelError,channelError,partitionError,partitionError,...impulseError,impulseError);save(fullfile(outDir,block_fir_state_matlab.mat),results,b);fprintf(BLOCK keep%.16g; conv%.16g; reset%.16g; state%.16g\n,...keepError,convError,resetError,stateError);fprintf(BLOCK tail%.16g; channels%.16g; partitions%.16g; impulse%.16g\n,...tailError,channelError,partitionError,impulseError);fprintf(BLOCK MATLAB PASS; blocks%d; taps%d; group_delay_samples%g\n,...numel(starts),L,(L-1)/2);endfunctionbmake_lowpass(L,cutoffCyclesPerSample)% Windowed-sinc FIR written directly; no fir1, hamming, or sinc dependency.m(-(L-1)/2:(L-1)/2);hzeros(L,1);z(m0);h(z)2*cutoffCyclesPerSample;h(~z)sin(2*pi*cutoffCyclesPerSample*m(~z))./(pi*m(~z));w0.54-0.46*cos(2*pi*(0:L-1)/(L-1));hh.*w;b(h/sum(h)).;endfunction[y,z,starts]run_blocks(x,b,pattern,keepState)% x is N-by-C, including C1; do not flatten a multi-channel matrix.yzeros(size(x));zzeros(numel(b)-1,size(x,2));startszeros(0,1);pos1;k1;whilepossize(x,1)countpattern(mod(k-1,numel(pattern))1);lastmin(size(x,1),poscount-1);if~keepState,z(:)0;end[part,z]filter(b,1,x(pos:last,:),z,1);y(pos:last,:)part;starts(end1,1)pos;%#okAGROWposlast1;kk1;endend7. 接到实际采集链路时先固定这几条规则初始化属于“流”不要放进每次 DMA 回调。每个通道有自己的状态状态要与采样流 ID 一起管理。只有接收完成、样本顺序确定后再推进滤波状态。重传包不能重复滤波丢包也不能假装未发生。缺样本以后补零、插值或重新启动都是明确的信号模型选择需要标记数据质量。系数改变、通道数改变或采样流重新开始时要定义状态如何迁移。直接沿用旧系数下的zf未必符合预期产品中常需要显式清零、过渡或交叉渐变。初次零状态仍然代表“记录之前的输入为零”。保留状态只消除了软件分块引入的重复启动不会恢复采集开始前未知的真实历史。与 FPGA 对齐时先对齐采样序号、16 点算法群延迟和 RTL 流水线延迟再讨论误差。硬件系数量化、舍入、饱和与握手停顿要单独建模。先用冲激、常量和几种不规则块长跑回归再上真实传感器波形。真实波形很难让人凭目视判断每个边界都正确。FIR 连续处理的验收标准很直接同一条输入流仅改变软件分块方式输出不应发生有意义的变化。把这个不变量写成自动测试比长期观察边界附近是否“看起来平滑”更可靠。延伸阅读MATLAB Hann 相干增益与补零标度对照。对时域状态连续性验收通过后再检查后续 FFT 的窗口与幅值标度避免把两个处理阶段的问题混在一起。工程问题如果你的 DMA 数据流中间丢了一个包你会选择补零后继续保留状态还是明确标记断流并重新初始化可以先用本文的不规则分块测试比较两种策略在丢包点附近的输出再按采集系统的数据质量要求决定。

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

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

免费获取报价 →
↑