Use Filter Constants to Hard Code Filter

75 次查看（过去 30 天）
Donald Hume
Donald Hume 2011-6-20
评论： Jan 2022-9-17
Hey,
I am trying to implement a real-time filter so am using MATLAB's butter() function to generate the needed [b,a] vectors
[b,a] = butter(4,4/35,'low');
Just to be clear I have used these generated vectors with the filter(a,b,data) function to successfully filter my data and it looks quite desirable. But as I am in the end trying to create a real time filter, I am trying to code this into a for-loop (for testing purposes). My code is as follows:
for n=5:1:length(x)
y(n) = b(1)*x(n)+b(2)*x(n-1)+b(3)*x(n-2)+b(4)*x(n-3)+b(5)*x(n-4)-a(2)*y(n-1)-a(3)*y(n-2)+a(4)*y(n-3)+a(5)*y(n-4);
end
This is the mathematical representation as far as I can gather from the doc: http://www.mathworks.com/help/techdoc/ref/filter.html
Can anyone tell me how I am incorrectly modeling the filter() command? I have also switched the a, b, column vectors (in case that was an issue). The above method just goes to infinity, and with a<->b the data just seems to be amplified.
Thanks for the help in advance.

请先登录，再回答此问题。

采纳的回答

Jan
Jan 2011-6-20
编辑：Jan 2014-10-26
The difference equation looks ok, but you do not show how e.g. "y(n-4)" is initialized.
Matlab's FILTER uses the "direct form II transposed" implementation, which is more efficient. Together with inital and final conditions:
function [Y, z] = myFilter(b, a, X, z)
% Author: Jan Simon, Heidelberg, (C) 2011
n = length(a);
z(n) = 0; % Creates zeros if input z is omitted
b = b / a(1); % [Edited, Jan, 26-Oct-2014, normalize parameters]
a = a / a(1);
Y = zeros(size(X));
for m = 1:length(Y)
Y(m) = b(1) * X(m) + z(1);
for i = 2:n
z(i - 1) = b(i) * X(m) + z(i) - a(i) * Y(m);
end
end
z = z(1:n - 1);
[EDITED]: A C-Mex implementation which handles arrays also: FEX: FilterM.
  19 个评论
显示 17更早的评论
Martin Goubej
Martin Goubej 2021-7-1
编辑：Martin Goubej 2021-7-1
@Francely Guzmán Otagri The provided code (as well as the Matlab filter command) implements the Transposed-Direct-Form-II structure of the IIR filter, check e.g. this page for explanation https://ccrma.stanford.edu/~jos/filters/Transposed_Direct_Forms.html#17507. The first line implements the upper summation term from the Fig. 9.2 and the rest deals with the update of the delay blocks.
As for the 'z' question - this comes form the filter structure, the uppermost summation term depends directly only on the output of the first internal delay. This is why it is assigned to the actual output value Y(m).
Jan
Jan 2022-6-8
With some more features:
function [Y, z] = myFilter(b, a, X, z)
% Author: Jan, Heidelberg, (C) 2022
na = numel(a); % Equilize size of a and b
n = numel(b);
if na > n
b(na) = 0;
n = na;
elseif na < n
a(n) = 0;
end
z(n) = 0; % Creates zeros if input z is omitted
b = b / a(1); % Normalize a and b
a = a / a(1);
Y = zeros(size(X));
if na > 1 % IIR filter
for m = 1:length(Y)
Y(m) = b(1) * X(m) + z(1);
for i = 2:n-1
z(i - 1) = b(i) * X(m) + z(i) - a(i) * Y(m);
end
z(n - 1) = b(n) * X(m) - a(n) * Y(m); % Omit z(n), which is 0
end
else % FIR filter: a(2:n) = 0
for m = 1:length(Y)
Y(m) = b(1) * X(m) + z(1);
for i = 2:n-1
z(i - 1) = b(i) * X(m) + z(i);
end
z(n - 1) = b(n) * X(m); % Omit z(n), which is 0
end
end
z = z(1:n - 1);
end

请先登录，再进行评论。

更多回答（3 个）

khatereh
khatereh 2012-1-6
Hi, I want to use your function instead of matlab filter function. I calculated the filter coefficient in b matrix and it is FIR filter so all the a values are 1. What should be the value for z in my case? I am confused how should I use z.
Thanks so much. Regards, KHatereh
  1 个评论
Jan
Jan 2014-10-26
编辑：Jan 2014-10-26
The meaning of z is explained in the documentation: doc filter.
The initial conditions for the internal state of the filter can be set such, that the transient effects are damped. Look into the code of the filtfilt function for a method to do this automatically.
Set z to zero, if you do not have any information about the signal.
For me the meaning of z got clear, when I examined this example: Imagine a long signal X, which is divided in 2 parts X1 and X2. Now the complete signal X is filtered with certain parameters and the initial settings z=0 (this means zeros(1,n-1) with n is the length of the filter parameters):
z = zeros(1, length(b) - 1);
Y = filter(b, a, X, z);
Now we do this for the first part:
z = zeros(1, length(b) - 1);
[Y1, z1] = filter(b, a, X1, z);
Now the output z1 is the internal state of the filter, a kind of history over the last elements. If we use the output z1 of the 1st part as input of the 2nd, we get exactly the same outpt as for the full signal:
Y2 = filter(b, a, X2, z1);
isequal(Y, [Y1, Y2]) % TRUE
But if we omit z1 as input for filtering X2, there is a small difference mostly at the start of Y2 due to the transient effects.
In this case, we do have some knowledge about the history of the internal filter state for X2, but for X1 this state is not defined and zeros are a fair guess, but not necessarily smart.

请先登录，再进行评论。

Yves
Yves 2018-5-10
Can someone please comment on whether this z/z1 or zi/zf - initial/final condition (delay) of the digital filter is equivalent to the state variables in state-space model (ABCD matrix) of the filter?
Michal
Michal 2022-2-21
Probably fully functional naive version of filter:
function [Y, z] = myFilter(b, a, X, z)
% a and b should have same order
na = length(a);
nb = length(b);
if na > nb
b = [b,zeros(1,na-nb)];
n = na;
elseif na < nb
a = [a,zeros(1,nb-na)];
n = nb;
else
n = na;
end
% naive filter implementation
z(n-1) = 0; % Creates zeros if input z is omitted
b = b / a(1); % normalize parameters
a = a / a(1);
Y = zeros(size(X));
for m = 1:length(Y)
Xm = X(m);
Y(m) = b(1) * Xm + z(1);
Ym = Y(m);
for i = 2:n-1
z(i - 1) = b(i) * Xm + z(i) - a(i) * Ym;
end
z(n - 1) = b(n) * Xm - a(n) * Ym;
end
end
  3 个评论
显示 1更早的评论
Michal
Michal 2022-6-8
编辑：Michal 2022-6-8
Your equilization length (a and b) method is good! But, the proposed modification:
z(n) = 0
broke the compatibility with original Matlab built-in function filter in a case of use z as an initial condition.
Error using filter
Initial conditions must be a vector of length max(length(a),length(b))-1,
or an array with the leading dimension of size max(length(a),length(b))-1
and with remaining dimensions matching those of x.
Jan
Jan 2022-9-17
Ran in:
@Michal: No, z(n)=0 does not break the compatibility with Matlab's filter, if it is applied inside the function and cropped finally: z = z(1:n - 1) as in my implementation above. The error message appears, if you expand z before calling filter, but this was not suggested.
Your method to expand z by the line:
z(n-1) = 0;
is a bug, because it sets the last element of z to 0. If z was [1,2,3] initially (then a and b have 4 elements), your code sets z to [1,2,0]. Simply try it:
x = rand(1, 100);
B = [0.000416546, 0.00124964, 0.00124964, 0.000416546];
A = [1, -2.68616, 2.41966, -0.730165]; % Lowpass Butterworth
z = [1,2,3];
y1 = filter(B, A, x, z);
y2 = michaelsFilter(B, A, x, z);
plot(y1, 'r')
hold('on')
plot(y2, 'b') % Obviously different
function [Y, z] = michaelsFilter(b, a, X, z)
% a and b should have same order
na = length(a);
nb = length(b);
if na > nb
b = [b,zeros(1,na-nb)];
n = na;
elseif na < nb
a = [a,zeros(1,nb-na)];
n = nb;
else
n = na;
end
% naive filter implementation
z(n - 1) = 0; % BUG! Must be z(n)=0, then crop z(n) after the loop.
% Or:
% if length(z) < n-1
% z(n - 1) = 0;
% end
b = b / a(1); % normalize parameters
a = a / a(1);
Y = zeros(size(X));
for m = 1:length(Y)
Xm = X(m);
Y(m) = b(1) * Xm + z(1);
Ym = Y(m);
for i = 2:n-1
z(i - 1) = b(i) * Xm + z(i) - a(i) * Ym;
end
z(n - 1) = b(n) * Xm - a(n) * Ym;
end
end
Do you see it now?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by