Re: filtering out harmonic trend

[Francesco Potorti`]

Francesco Potortì:
>> I have a 200,000 samples measurement, which is a sort of noise with
>> bandwidth in the order of the inverse of 500 samples.  Superimposed to
>> this noise there is a very slow sinusoid with a period of about 300,000
>> samples.  This means that my measurements do not even see a complete
>> cycle of this sinusoid.
>>
>> I want to estimate the sinusoid (in order to remove it from the
>> measurement).  So I need some sort of very-low pass filter.  However,

Dmitri A. Sergatskov:
>If you know the frequency of the sinusoid (e.g. if it is a power line
>interferance), then you should be able to do a generalized linear
>fit to :
>
>Y = A*sin(Wt) + B*cos(Wt) + C

Yes, that was my first idea.  But using leasqr with 200,000 samples I
fear it takes a lot of time.  I have not tried this yet.  Thank you.

Rafael Laboissiere:
>filtfilt

After looking at the list of functions in octave-forge, this was my
second idea, but since I had never heard of it, I was doubtful.  Thank
you.

I tried filtfilt, and found that not being able to specify initial and
final states is a big drawback for my application.  So I patched
filtfilt to add initial and final states.  It works, but it can be
improved.  There are three issues:

1) it adds two optional parameters, which are not present in Matlab's
   version
2) doing like Matlab would be more compatible and easier to use
3) it gives results different from the current implementation.

About the second issue, I see that, in the source, reference is made to
filtic.  I tried, but because of little time I could not understand how
does it work.  I figure that it would be very easy to change the values
of si and sf in the code below to use filtic.  If only someone could
give me a hint on how to use it, I would gladly improve on the patch.

--- /usr/share/octave/site/api-v13/m/octave2.1-forge/signal/filtfilt.m~ 
2006-11-28 19:28:15.000000000 +0100
+++ /usr/share/octave/site/api-v13/m/octave2.1-forge/signal/filtfilt.m  
2007-12-04 18:47:17.000000000 +0100
@@ -21,13 +21,7 @@
 ## magnitude response in the process. That's the theory at least.  In
 ## practice the phase correction is not perfect, and magnitude response
 ## is distorted, particularly in the stop band.
-##
-## In this version, I zero-pad the end of the signal to give the reverse
-## filter time to ramp up to the level at the end of the signal.
-## Unfortunately, the degree of padding required is dependent on the
-## nature of the filter and not just its order, so this function needs
-## some work yet.
-##
+####
 ## Example
 ##    [b, a]=butter(3, 0.1);                   % 10 Hz low-pass filter
 ##    t = 0:0.01:1.0;                         % 1 second sample
@@ -56,18 +50,24 @@
 ## TODO:   Don't know for what n this would be faster, if any, but the
 ## TODO:   memory saving might be nice.
 
-function y = filtfilt(b, a, x)
-  if (nargin != 3)
-    usage("y=filtfilt(b,a,x)");
+
+function y = filtfilt(b, a, x, si, sf)
+  if (nargin < 3 || nargin > 5)
+    usage("y=filtfilt(b,a,x,si,sf)");
   end
 
+  slen=max(length(a),length(b)); # size of status
+  if (nargin < 4)               # no initial state provided
+    si = zeros(1,slen);
+  endif
+  if (nargin < 5)               # no final state provided
+    sf = zeros(1,slen);
+  endif
   if (rows(x) == 1)
-    y = filter(b,a,[x, zeros(1,2*max(length(a),length(b)))]);
-    y = fliplr(filter(b,a,fliplr(y))); 
-    y = y(1:length(x));
+    y = filter(b,a,x,si);
+    y = fliplr(filter(b,a,fliplr(y),sf)); 
   else
-    y = filter(b,a,[x ; zeros(2*max(length(a),length(b)), columns(x))]);
-    y = flipud(filter(b,a,flipud(y))); 
-    y = y(1:rows(x),:);
+    y = filter(b,a,x,si);
+    y = flipud(filter(b,a,flipud(y),sf)); 
   endif
 endfunction

-- 
Francesco Potortì (ricercatore)        Voice: +39 050 315 3058 (op.2111)
ISTI - Area della ricerca CNR          Fax:   +39 050 315 2040
via G. Moruzzi 1, I-56124 Pisa         Email: address@hidden
Web: http://fly.isti.cnr.it/           Key:   fly.isti.cnr.it/public.key