MATLAB Answers

Johannes
0

Filter with variable span

Asked by Johannes
on 29 May 2018
Latest activity Edited by Anton Semechko on 31 May 2018
Hello everybody,
I am in need of a filter with a variable span as my data points (x-values) have non-uniform spacing. I would like to use a moving average filter and/or Savitzky Golay Filter. However, the problem is that the spacing of my x-values is not uniform but decreasing (at least it is monotonously decreasing, but their is no systematic trend). When I use a simple MA or SG filter, I do have the problem that the data in the regions of low point density (large spacing, in the end) is smoothed way too much compared to the data in the beginning (high point density). Thus, I would need a filter function which automatically adjusts to that (diff(x) so to say). A least for the MA filter, I could think of a self-made solution using simple loops. But I guess that's very inefficient and I might not be the first one to seek for such a solution.
Any help is therefore much appreciated.
Thank you very much, Johannes

  0 Comments

Sign in to comment.

3 Answers

Answer by Anton Semechko on 29 May 2018

Simplest solution would be to resample your signal on a uniformly-spaced grid using linear interpolation (see 'interp1' function) and then feed resampled signal to the filter of your choice

  0 Comments

Sign in to comment.


Answer by Johannes
on 29 May 2018

Hi Anton, thank you for your answer. I also thought about the same thing. The grid would then have to be as fine as the highest point density I have in my data. However, I was afraid that this might mess with the intrinsic white noise in my data. Or wouldn't you expect this to be a problem? Thanks!

  1 Comment

If you use linear interpolation, it will not introduce any spurious noise into your signal. There are other options, however, such as Laplacian or bilateral filtering , that operate directly on non-uniformly sampled signals. The main idea there is to use filter weights that are inversely related to the distances among the nodes.

Sign in to comment.


Answer by Anton Semechko on 29 May 2018
Edited by Anton Semechko on 31 May 2018

Here is an example of a Laplacian filtering function I just put together:
N=5E3; % number of samples
X=sort(rand(N,1)); % sampling locations
Fo=sin(40*X)+cos(10*X); % signal
F=Fo+0.75*(2*rand(N,1)-1); % signal + white noise
% Filter the signal
c=[0.1 0.9]; % filter parameters; c(1) and c(2) penalize derivatives of orders greater than 1 and 2, resp
t=1E5; % diffusion time; greater values produce smoother signal
vis=true; % visualize output
[Fout,H]=LaplacianFilt_1D(X,F,c,t,vis);
% Show the original signal
h3=plot(X,Fo,'--g','LineWidth',2);
h=legend([H h3],{'input' 'output' 'original'});
set(h,'FontSize',20,'Location','NorthEast')
Main function:
[Fout,H]=LaplacianFilt_1D(X,F,c,t,vis)
% Smooth a non-periodic 1D signal with non-uniform sampling intervals using
% Laplacian (i.e., diffusion) filtering.
%
% - X : N-by-1 vector of uniformly increasing sample "locations"
% - F : N-by-1 vector of sample values corresponding to X
% - c : c=[c1 c2], where c1 and c2 are Laplacian and bi-Laplacian
% weights
% - t : time step
% - vis : set vis=true to visualize output
%
% AUTHOR: Anton Semechko (a.semechko@gmail.com)
%
if nargin<2 || isempty(X) || isempty(F)
error('Insufficient number of input arguments')
end
X=X(:);
F=F(:);
N=numel(X);
if N~=numel(F)
error('Invalid entry for 2nd input argument (F)')
end
if nargin<3 || isempty(c)
c=[0.1 0.9];
elseif numel(c)~=2 || ~isnumeric(c) || sum(c<0)
error('Invalid entry for 3rd input argument (c)')
end
c=c+eps;
c=c/sum(c);
if nargin<4 || isempty(t)
t=10;
elseif numel(t)~=1 || ~isnumeric(t) || t<eps
error('Invalid entry for 4th input argument (t)')
end
if nargin<5 || isempty(vis)
vis=false;
elseif numel(vis)~=1 || ~logical(vis)
error('Invalid entry for 5th input argument (vis)')
end
% Fujiwara Laplacian
W=X(2:N)-X(1:(N-1)); % distance between the nodes
W_min=max(min(W)/1E3,eps);
W=max(W,W_min);
Ew=max(W)./W; % edge weights
E=[1:(N-1);2:N]';
Ne=N-1;
i=repmat((1:Ne)',[2 1]);
j=E(:);
s=ones(2*Ne,1); s((Ne+1):end)=-1;
A=sparse(i,j,s,Ne,N,2*Ne);
L=A'*spdiags(Ew,0,Ne,Ne)*A;
L=(N/trace(L))*L;
% Filter
A=speye(N)+t*(c(1)*L+c(2)*L^2);
Fout=A\F;
% Visualize
H=[];
if ~vis, return; end
figure('color','w')
h1=plot(X,F,'-k','LineWidth',2);
hold on
h2=plot(X,Fout,'-r','LineWidth',3);
H=[h1 h2];

  0 Comments

Sign in to comment.