Loop in filtfilt not really necessary ?
Deric Wisleder
is the nested for loop on line 48
line 46-53 of the filtfilt function:
if (n>1)&(m>1)
y = x;
for i=1:n % loop over columns
y(:,i) = filtfilt(b,a,x(:,i));
end
return
% error('Only works for vector input.')
end
A colleague replaced the reflected extrapolations in filtfilt with
spline extrapolations as follows:
% Extrapolate beginning and end of data sequence using a "seventh
order polynomial
% method". Slopes of original and extrapolated sequences match at
% the end points.
% This reduces end effects.
ext2=8 ; ext=20 ;
q=(0:ext2-1) ;
q2=(-1*(ext2)+1:1:0) ;
p=polyfit(q,x(1:ext2)',7) ;
p2=polyfit(q2,x(len-ext2+1:len)',7) ;
q3=(-1*ext:1:-1)/50 ; q4=(1:ext)/50 ;
d=polyval(p,q3) ; d2=polyval(p2,q4) ;
y = [d';x;d2'];
polyval reportedly works on matrices, but I do not know about
polyfit. I can use vertcat to vectorize appending of the
extrapolated segments. Thus, polyfit and polyval seem to be the only
commands that need to be vectorized to remove the loop in my current
version of filtfilt. Likewise the linear extrapolation line in the
standard version of filtfilt would appear to be the only line to
vectorize for removal of the loop:
y =
[2*x(1)-x((nfact+1):-1:2);x;2*x(len)-x((len-1):-1:len-nfact)];
“Vectorize “filtfilt” recursion” got no response. Does this post
clarify the question and isolate the problem? Can someone please
help me to remove the loop in filtfilt if possible.
Thank you very much.
Deric Wisleder
Matthew Carson
was very slow. Keep in mind that my data all begins and ends with
zero, and I did a simple moving average filter with
a=1;
b=ones(1,(num_avg+1)/2)/((num_avg+1)/2);
I was able to replace
filtfilt(b,a,signal)
with
filt_sig=filter(b,a,signal);
filt_sig=filt_sig(720:-1:1,:,:);
filt_sig=filter(b,a,filt_sig);
filt_sig=filt_sig(720:-1:1,:,:);
I got exactly the same results as filtfilt and it took 0.05 seconds
instead of 0.300 seconds.
D Wisleder
>
> I was having a problem where calling filtfilt with a 720x30x6
>
> I was able to replace
> filtfilt(b,a,signal)
>
>
> I got exactly the same results as filtfilt
> and it took 0.05 seconds instead of 0.300
> seconds.
Thank you Matthew,
I use a butterworth filter in accordance with other researchers in my
area. If all matrix columns are to be filtered with the same cutoff,
the filter coefficients are calculated in three simple lines before
calling filtfilt (my developing version). Your replacement lines are
similar to the sequence in filtfilt (filter, reverse data, filter
again, and reverse data again) and could be implemented there. It
still seems that the only purpose for the loop lies in creating and
removing extrapolated segments.
Can someone help with vectorizing polyfit (Does it already work on
independent columns?) and polyval, andthe one-line padding in the
standard filtfilt? Thank you.
Deric
Matthew Carson
tell. Actually polyval will work on a matrix but seems to not be
able to use a matrix of coefficients. It looks like it only uses the
first column.
iow
x=rand(10,10);
a=rand(10,10);
b=polyval(a,x);
c=polyval(a(:,1),x);
b==c will be all 1's
you can always do a loop just on the polyfit:
ext2=8 ; ext=20 ;
q=(0:ext2-1) ;
q2=(-1*(ext2)+1:1:0) ;
q3=(-1*ext:1:-1)/50 ; q4=(1:ext)/50 ;
for n=1:length(data)
p=polyfit(q,x(1:ext2)',7) ;
p2=polyfit(q2,x(len-ext2+1:len)',7) ;
d(n,:)=polyval(p,q3);
d2(n,:)=polyval(p2,q4) ;
end
y = [d';x;d2'];
then filter the whole thing together
D Wisleder
> >
> polyfit and polyval do not seem to work on matrices
That is too bad.
> you can always do a loop just on the polyfit:
> then filter the whole thing together
Your suggestion was very helpful. My implementation turned out a
little different, but this portion seems to work:
ext2=8 ; ext=20 ;
q=(0:ext2-1) ;
q2=(-1*(ext2)+1:1:0) ;
q3=(-1*ext:1:-1)/50 ;
q4=(1:ext)/50 ;
ldr = ones(20,length(x(1,:))) ; trlr = ldr ; % Pre-allocate
for cols2=1:length(x(1,:)) ;
p=polyfit(q,x(1:ext2,cols2)',7) ;
p2=polyfit(q2,x(len-ext2+1:len,cols2)',7) ;
ldr(:,cols2)=polyval(p,q3)' ;
trlr(:,cols2)=polyval(p2,q4)' ;
end
y = [ldr; x; trlr] ; % Vertical concatination
Do I need to adjust y(1) to y(1,:) in the filter commands below?
% filter, reverse data, filter again, and reverse data again
y = filter(b,a,y,[zi*y(1)]) ;
y = flipud(y) ;
y = filter(b,a,y,[zi*y(1)]) ;
y = flipud(y) ;
Matthew Carson
depending on the direction of the vectors. It should end up a
matrix. The help says
Input zi is 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.
It only talks about using zi in conjunction with the zf output
however, and I see you don't have that...
Also I tried flipud versus the indexing I put in there and found the
indexing to be slightly faster. YMMV
> % filter, reverse data, filter again, and reverse data again
> y = filter(b,a,y,[zi*y(1)]) ;
> y = flipud(y) ;
> y = filter(b,a,y,[zi*y(1)]) ;
> y = flipud(y) ;
-Matt
Deric Wisleder
Thank you very much! I have this working and can post or attach
materials if desired.
Matthew Carson wrote:
>
> I think you would need to have zi*y(1,:), or possibly zi'*y(1,:)
> depending on the direction of the vectors. It should end up a
> matrix.
> The help says
> Input zi is 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.
>
> It only talks about using zi in conjunction with the zf output
> however, and I see you don't have that...
My data matrices are column vectors, so I indexed initial conditions
accross the columns but left the zi transform zi*y(1,:) as in the
original filtfilt.
> Also I tried flipud versus the indexing I put in there and found
> the
> indexing to be slightly faster. YMMV
>> % filter, reverse data, filter again, and reverse data again
>> y = filter(b,a,y,[zi*y(1)]) ;
>> y = flipud(y) ;
>> y = filter(b,a,y,[zi*y(1)]) ;
>> y = flipud(y) ;
I switched back to your indexing; although, I have not used the
profiler yet (on any code).
Thank you again