我如何自己编写Matlab的“过滤器"功能? [英] How can I write the Matlab "filter"-function myself?

问题描述

我想在1D信号上使用Butterworth滤波器.在Matlab中，脚本如下所示:

I would like to use a Butterworth filter on a 1D-Signal. In Matlab the script would look like this:

 f=100;
 f_cutoff = 20; 
 fnorm =f_cutoff/(f/2);
 [b,a] = butter(8,fnorm,'low');
 filteredData = filter(b,a,rawData); % I want to write this myself

现在，我不想直接使用Matlab中提供的 filter 函数，而是自己编写. 在Matlab文档中，其描述如下:

Now I don't want to directly use the filter-function given in Matlab but write it myself. In the Matlab documentation it's described as follows:

过滤器功能实现为直接形式II换位结构，

The filter function is implemented as a direct form II transposed structure,

y(n)= b(1)* x(n)+ b(2)* x(n-1)+ ... + b(nb + 1)* x(n-nb) -a(2)* y(n-1)-...-a(na + 1)* y(n-na)

y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb) - a(2)*y(n-1) - ... - a(na+1)*y(n-na)

其中n-1是滤波器阶数，可同时处理FIR和IIR滤波器[1]，na是反馈滤波器阶数，nb是前馈滤波器阶数.

where n-1 is the filter order, which handles both FIR and IIR filters [1], na is the feedback filter order, and nb is the feedforward filter order.

所以我已经尝试过编写这样的函数:

So I've already tried to write the function like that:

f=100;
f_cutoff = 20; 
fnorm =f_cutoff/(f/2);
[b,a] = butter(8,fnorm,'low');
for n = 9:size(rawData,1)
    filteredData(n,1) = b(1)*n + b(2)*(n-1) + b(3)*(n-2) + b(4)*(n-3) + b(5)*(n-4) ...
                      - a(2)*rawData(n-1,1) - a(3)*rawData(n-2,1) - a(4)*rawData(n-3,1) - a(5)*accel(n-4,1);
end

但这不起作用.你能帮我么?我在做什么错了?

But that's not working. Can you please help me? What am I doing wrong?

此致， 塞尔多(Cerdo)

Sincerely, Cerdo

PS:可以在此处找到过滤器文档:展开更多关于->算法

PS: the filter documentation can be foud here: http://www.mathworks.de/de/help/matlab/ref/filter.html#f83-1015962 when expanding More About -> Algorithms

推荐答案

我找到了描述在Matlab过滤器函数中使用的Direct Form II Transposed的文本，并且效果很好.请参见下面的脚本.其他实现也可用，但错误大约为1e-15，您可以通过自己运行脚本来查看.

I have found a text described the Direct Form II Transposed used in the Matlab filter function and it works perfectly. See script below. Other implementations are also available but with error of around 1e-15, you'll see this by running the script yourself.

%% Specification of the Linear Chebysev filters
clc;clear all;close all
ord = 5; %System order (from 1 to 5)
[bq,aq] = cheby1(ord,2,0.2);theta = [bq aq(2:end)]';
figure;zplane(bq,aq); % Z-Pole/Zeros
u = [ones(40,1); zeros(40,1)];
%% Naive implementation of the basic algorithm
y0 = filter(bq,aq,u); % Built-in filter
b = fliplr(bq);a = fliplr(aq);a(end) = [];
y1 = zeros(40,1);pad = zeros (ord,1);
yp = [pad; y1(:)];up = [pad; u(:)];
for i = 1:length(u)
    yp(i+ord) = sum(b(:).*up(i:i+ord))-sum(a(:).*yp(i:i+ord-1));
end
y1 = yp(ord+1:end); % Naive implementation
err = y0(:)-y1(:);
figure
plot(y0,'r')
hold on
plot(y1,'*g')
xlabel('Time')
ylabel('Response')
legend('My code','Built-in filter')
figure
plot(err)
xlabel('Time')
ylabel('Error')
%% Direct Form II Transposed  
% Direct realization of rational transfer functions
% trps: 0 for direct realization, 1 for transposed realisation 
% b,a: Numerator and denominator
% x:   Input sequence
% y:   Output sequence
% u:   Internal states buffer

trps = 1;
b=theta(1:ord+1);
a=theta(ord+2:end);
y2=zeros(size(u));
x=zeros(ord,1);
%%
if trps==1
    for i=1:length(u)
        y2(i)=b(1)*u(i)+x(1);
        x=[x(2:ord);0];
        x=x+b(2:end)*u(i)-a*y2(i);
    end
else
    for i=1:length(u)
        xnew=u(i)-sum(x(1:ord).*a);
        x=[xnew,x];
        y2(i)=sum(x(1:ord+1).*b);
        x=x(1:ord);
    end
end
%%
err = y2 - filter(bq,aq,u);
figure
plot(y0,'r')
hold on
plot(y2,'*g')
xlabel('Time')
ylabel('Response')
legend('Form II Transposed','Built-in filter')
figure
plot(err)
xlabel('Time')
ylabel('Error')
% end

这篇关于我如何自己编写Matlab的“过滤器"功能?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！