FiltFilt function initial and final conditions

46 次查看（过去 30 天）
Brecht Vermeulen
Brecht Vermeulen 2018-8-29
编辑： Bruno Luong 2018-10-16
Hi,
I would like to filter a signal using filtfilt in order to have:
- zero phase - to have the initial and final sample of the filtered data the same as the raw data
Following this link: http://www.mechanicalvibration.com/filtfilt_Casual_versus_non_.html#foot60 the initial and final point are maintained using filtfilt. However when I do this, they are not.
The commands I use are: b = [0.0675, 0.1349, 0.0675]; % Low pass Butterworth filter a = [1.0, -1.1430, 0.4128]; Sig_filtered=filtfilt(b, a, Sig_raw);
Can anyone explain this to me or point me towards a filter function that can do this?
Thank you in advance.
Kr, Brecht

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

采纳的回答

Jan
Jan 2018-8-29
As written in the documentation, and as you can see inside the code, filtfilt mirrors the initial and final values to reduce the transitional effects. This cannot preserve the marginal values in every case.
It is not clear, what you exactly want to achieve. Please explain uniquely and clearly what you need.
  6 个评论
显示 4更早的评论
Jan
Jan 2018-10-16
编辑：Jan 2018-10-16
@Robert: In the documentation of filtfilt.m you find the source of the applied method:
Frederik Gustafsson, Determining the initial state in forward-backward filtering, IEEE Transactions on Signal Porcessing, pp 998-992, April 1996, Volume 44, Issue 4
Therefore I disagree with the opinion, that MathWorks is sloppy with the exact description of the mathematical decisions.
From a physical point of view, the filter is a resonator. The filtered output depends on the input signal and the state of the internal "oscillators". Then the method of reducing the transitional effects in filtfilt can be interpreted such, that it is assumed, that the signal can be expanded by subtracting the flipped initial phase from the first element. This is a fair guess for a sine wave, for a static signal, for pure noise. The spectrum of the signal is conserved approximately and this does not create a DC offset. But of course e.g. for a square or triangle wave the results are not perfect. There cannot be a perfect method in general. As soon as you do have knowledge about the real signal in before and after the filtered section, it is wise to use it.
Is this "built in" or is there source to this function?
It seems like you have found out already, that the code is shipped with Matlab. Simply try
type filtfilt.m
You find the code there and the above mentioned reference for the algorithm. This reference is mentioned in the documentation also: doc filtfilt (link)
Robert Bristow-Johnson
Robert Bristow-Johnson 2018-10-16
编辑：Robert Bristow-Johnson 2018-10-16
Bruno, TMW and Cleve know that i am an historic pain-in-arse regarding MATLAB and will continue to be a grumpy old man as long as the indexing origin is hard-wired to 1 (i still can't fathom that, such a basic and fundamental problem). but now i have at least a couple of audio friends that work for TMW and they gave me a cute little T shirt so now i am less critical than i was on comp.dsp and comp.soft-sys.matlab in the '90s and early 2000s.
but even today, it's disgusting. my MATLAB code is disgusting. ugliest friggin' code i have ever written. it's all ugly ugly ugly. and that 1 origin convention is the main reason. not being able to have negative indices is another reason. and the order of coefficients for polyval() and polyfit() is wrong.
so i guess i am still a grumpy old DSPer.

请先登录，再进行评论。

更多回答（3 个）

Bruno Luong
Bruno Luong 2018-10-16
编辑：Bruno Luong 2018-10-16
Here is a demo of using filtfilt on periodic data
  • The first method use filtfilt on one period
  • The second method stitches rawdata one period on the left and one on the right, call filtfilt on the three periods, then crop the middle period (Jan's suggestion)
  • The third method does like the second but smoothly interpolate in a transition of 0.1 period.
As you can see the first one have a problem of stitching at the boundary, the second and third methods give almost the same result, but I know the third one must have a smooth transition.
% True Periodic signal
x = linspace(0,2*pi,100);
y = sin(x) - 0.5*sin(2*x+1) + 0.6*sin(3*x+2) - 0.3*sin(4*x+4);
% Add noise
sigma = 0.3;
y = y + sigma*randn(size(y));
% Filter design
[b,a] = butter(8,0.1);
% First method
yff = filtfilt(b,a,y);
% Second & Third methods
nperiods = 3;
[xd,yd] = duplicate(x,y,nperiods);
y123 = filtfilt(b,a,yd);
y1 = crop(y123,nperiods,1);
y2 = crop(y123,nperiods,2);
% Second method
yt = y2;
% Third method
w = max(0,interp1(x(end)*[0.9 1],[0,1],x,'linear','extrap'));
yi = (1-w).*y2 + w.*y1;
% Generate graphics
close all
hold on
xp = x/(2*pi);
h(1)=plotdup(xp,y,'color',0.7*[1 1 1]);
h(2)=plotdup(xp,yff,'b');
h(3)=plotdup(xp,yt,'g');
h(4)=plotdup(xp,yi,'r');
xline(1);
xline(2);
legend(h,'noisy data','fiiltfilt','ff-truncated','ff-interpolate')
function [xd,yd] = duplicate(x,y,nperiods)
if nargin < 3
nperiods = 3;
end
x = x(:);
dx = x(end)-x(1);
xd = x(1:end-1) + dx*(0:nperiods-1);
xd = [xd(:)', dx*(nperiods-1)+x(end)];
yd = [repmat(y(1:end-1),1,nperiods), y(end)];
end
function y = crop(yd,nperiods,i)
yend = yd(end);
yd = reshape(yd(1:end-1),[],nperiods);
if i < nperiods
yend = yd(1,i+1);
end
y = [yd(:,i); yend]';
end
function h = plotdup(x,y,varargin)
[xd,yd] = duplicate(x,y);
h = plot(xd,yd,varargin{:});
end
Bruno Luong
Bruno Luong 2018-10-16
编辑：Bruno Luong 2018-10-16
Now a rigorous approach of zero-phase FIR on circular data. One need to build a circular matrix from the coefficients A and B, the similar to filt-filt apply the circular-filter in both directions to undo the phase-shift.
I call it fourth method, a ff-circular. The code is following
% True Periodic signal
x = linspace(0,2*pi,100);
y = sin(x) - 0.5*sin(2*x+1) + 0.6*sin(3*x+2) - 0.3*sin(4*x+4);
% Add noise
sigma = 0.3;
y = y + sigma*randn(size(y));
% Filter design
[b,a] = butter(8,0.1);
% First method
yff = filtfilt(b,a,y);
% Fourth method
n = length(y)-1;
Sa = cmat(a,n);
Sb = cmat(b,n);
yc=Sa\(Sb*y(1:end-1)');
yc=Sa\(Sb*flip(yc));
yc = flip(yc([end 1:end]))';
% Generate graphics
close all
hold on
xp = x/(2*pi);
h(1)=plotdup(xp,y,'color',0.7*[1 1 1]);
h(2)=plotdup(xp,yff,'b');
h(3)=plotdup(xp,yc,'r');
xline(1);
xline(2);
legend(h,'noisy data','fiiltfilt','ff-circular')
function S = cmat(a, n)
A = repmat(flip(a),n);
d = 0:length(a)-1;
S = spdiags(A,d,n,n) + ...
spdiags(A,d-n,n,n);
end
function [xd,yd] = duplicate(x,y,nperiods)
if nargin < 3
nperiods = 3;
end
x = x(:);
dx = x(end)-x(1);
xd = x(1:end-1) + dx*(0:nperiods-1);
xd = [xd(:)', dx*(nperiods-1)+x(end)];
yd = [repmat(y(1:end-1),1,nperiods), y(end)];
end
function h = plotdup(x,y,varargin)
[xd,yd] = duplicate(x,y);
h = plot(xd,yd,varargin{:});
end
Brecht Vermeulen
Brecht Vermeulen 2018-9-4
I expected filtfilt to maintain exactly the initial and final sample. When heavily filtering this deviation might become quite big.
What I exactly needed is filtering of a cyclic signal, where the start and endpoint were exactly the same. I manually have written an adjustment code that did this.
Thank you for the clarification.
  4 个评论
显示 2更早的评论
Bruno Luong
Bruno Luong 2018-10-16
编辑：Bruno Luong 2018-10-16
@Jan: A better way than crude cropping is interpolating the filtered signal over overlapping zone (imagine the signal wrapped on the circle) by a weighted of two signals that make a smooth transitions with continuity, the weight selected to can be infinity smooth if such smoothness is desired.
Jan
Jan 2018-10-16
编辑：Jan 2018-10-16
@Bruno: I agree. I'm only not satisfied with such photoshopping of measurement data. Smoothing a signal until it fits the expectations exactly is a construction of a result and not necessarily a measurement. Another problem could be a high frequency signal near to the Nyquist frequency. Then a minimal timeshift can cause a cancellation of the signal. So cropping is less smart, but perhaps this crudeness preserves the nature of the original signal.
As usual in signal processing, filtering is not trivial. You cannot expect that filter or filtfilt remove the noise magically and a deep knowledge of the effects of filtering is required. A pure noise can be filtered, until a clean periodic sine-wave is found as "result", but this has no scientific power anymore.
But maybe the OP is aware of such problems. It was only "heavy filtering" combined with "maintain exactly the initial and final sample" which let me mention the potential overdoing at filtering.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Smoothing and Denoising 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by