Fortran implementation of Matlab's filtfilt

Hi All,

I'm hunting around for a head start on implementing a zero phase
distortion filter (a.k.a "filtfilt" in Matlab) using Fortran.  I'm
familiar with Matlab but not so much with Fortran - I essentially need to
code something very similar to filtfilt using Fortran so that the
filtering routine can be called from other Fortran based code (using
static filter coefficients as designed in Matlab).

Any pointers appreciated!.

Regards,
Christopher.





_____________________________________
Do you know a company who employs DSP engineers?  
Is it already listed at http://dsprelated.com/employers.php ?

On 4 May, 13:35, "cwatson" <cwat...@utas.edu.au> wrote:
> Hi All,
>
> I'm hunting around for a head start on implementing a zero phase
> distortion filter (a.k.a "filtfilt" in Matlab) using Fortran.  I'm
> familiar with Matlab but not so much with Fortran - I essentially need to
> code something very similar to filtfilt using Fortran so that the
> filtering routine can be called from other Fortran based code (using
> static filter coefficients as designed in Matlab).

First of all, the 0-phase filter is non-causal, and can therefore only
be
used in off-line or batch mode. This means that you have all the data
you want to filter available when you start the filtering operation.

You can proceed in several ways:

One method is to

- Compute the DFT of the data
- Compute the DFT of the filter
- Multiply the data spectrum with the *magnitude* of the filter
spectrum,
  i.e. the phase = 0 for all coefficients in the filter spectrum.
- Compute the IFFT of the resulting spectrum.

Another method is to

- Apply the filter in the forward direction of the data
- Flip the data around and apply the same filter on the time-reversed
  data.

In the special case when the filter is a symmetric FIR filter with an
odd number of coefficients N, you can assign the output to index
(N-1)/2 instead of index N-1, which is the normal way.

Rune

On May 4, 6:35 am, "cwatson" <cwat...@utas.edu.au> wrote:
> Hi All,
>
> I'm hunting around for a head start on implementing a zero phase
> distortion filter (a.k.a "filtfilt" in Matlab) using Fortran.  I'm
> familiar with Matlab but not so much with Fortran - I essentially need to
> code something very similar to filtfilt using Fortran so that the
> filtering routine can be called from other Fortran based code (using
> static filter coefficients as designed in Matlab).
>
> Any pointers appreciated!.
>
> Regards,
> Christopher.
>


I think this is already done in many books, such as
"Numerical Recipes in Fortran".

Julius

On 2007-05-04, cwatson <cwatson@utas.edu.au> wrote:
> I'm hunting around for a head start on implementing a zero phase
> distortion filter (a.k.a "filtfilt" in Matlab) using Fortran.

The octave-forge version is just a filter with a pad chosen based on
the filter order followed by a filter in the opposite direction:

	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));

Maybe Matlab's is more clever wrt to overall gain?

-- 
Ben Jackson AD7GD
<ben@ben.com>
http://www.ben.com/