[Matlab] 滤波器filter函数造轮子及使用代码生成进行速度优化

virter 2021-08-02 原文

之前做脑机接口上位机的时候需要对数据进行实时滤波，也就是需要对数据进行分段滤波，保存滤波器前一次的历史状态。翻了翻MATLAB官方文档的filter函数发现有这个功能，不过他们的函数说明是用相位及延迟进行设置，看了半天没理解什么意思，为了保险起见我自己造了个轮子实现简单的IIR滤波器。filter函数的官方文档：https://ww2.mathworks.cn/help/matlab/ref/filter.html?s_tid=doc_ta 在官方文档中IIR滤波器的详细原理都已经讲的很明白了不再煞述。自己造的轮子也很简单，关于滤波器的参数都被保存在myIIR的结构体中，包括滤波器系数a,b，历史数据x,y，以及滤波后的数据FilteredData。在实时检测系统中一般对速度由较强的要求，所以我们可以将MATLAB代码转成mex文件进行提速。不是所有MATLAB代码都是可以转成mex或C的，总结来说如果要想顺利转成C的话还是得用C的思想写MATLAB代码，MATLAB那些奇技淫巧能少用就少用，要用某些函数的手先到官方文档上查一下有没有相应的C/C++ Coder的支持。这样这个简单的滤波函数代码如下：

function myIIR = myFilter(data,myIIR)
    for iSample = 1:size(data,2)
        myIIR.x = circshift(myIIR.x,1,2);
        myIIR.x(:,1) = data(:,iSample);
%         myIIR.x = [data(:,iSample),myIIR.x];    % PUSH
%         myIIR.x(:,end) = [];    % POP
        y = sum(myIIR.b.*myIIR.x,2)-sum(myIIR.a.*myIIR.y,2);
        myIIR.y = circshift(myIIR.y,1,2);
        myIIR.y(:,1) = y;
%         myIIR.y = [y,myIIR.y];  % PUSH
%         myIIR.y(:,end) = [];    % POP
        myIIR.FilteredData(:,iSample) = y;
    end
end

代码注释掉的部分就是代码生成不支持的语句，用循环位移函数circshift函数代替了。为了对比，内置函数也写成这种结构体的形式，代码如下：

function IIR = Filter(data,IIR)
    [IIR.FilteredData,IIR.zf] = filter(IIR.b,IIR.a,data,IIR.zf,2);
end

进行这样设置以后，可以进行mex文件生成，生成后在本机上进行测试。测试数据是一个cell元胞数组，包含10000个8×256矩阵，模拟10000次8通道滤波。本机测试环境：

程序运行的MATLAB版本为MATLAB-2018b，我分别在MATLAB-2016a, -2018a,-2018b三个不同版本的MATLAB上生成了3种mex文件对比不同版本MATLAB的差异，结果如下：

可以看出来自己用MATLAB造的filter函数的轮子远远没有内置函数的快，但转为mex文件以后可以提高至少10倍的速度。而内置函数一直是最快的，甚至生成mex文件后速度还有一定的下降。不同版本之间生成的代码也是有差异的，2016a中生成的轮子的mex明显比2018a和2018b中生成的要快，但内置函数转为mex后反而慢了。。这个现象也是挺奇葩

综上，我感觉除非真的对速度有较为苛刻的要求，否则还是用MATLAB语法里的奇技淫巧写代码效率比较高，若要转为mex的话需要参考不少C 的语法，那MATLAB代码写着就没啥意思了