FILTER code (FORTRAN)

Fred Marshall schrieb:
> I once was able to speed up a "filter" by better than a factor of 10 by 
> going to the frequency domain.  It's what was predicted and it's what we 
> got.
> 
> The key is that the filter and the data are relatively long sequences.  The 
> longer the sequences, the more computational saving.
> 
> Here's how it's done:
> 
> 1) Determine a FIR filter.  This is a step unto itself.  Call the length of 
> the filter "N".
> 
> 2) Determine the length of data sequence that you will process at any one 
> time.  Call the length of the data sequence "M".
> 
> 3) Calculate M + N -1.  This, or a greater value, will be the length of the 
> FFT / IFFT sequences.  Decide on this value "R" - maybe a power of 2 if your 
> FFT is finicky about that.
> 
> 3) Pad the filter with zeros to length "R".
> 
> 4) Pad the data with zeros to length "R".
> 
> 5) FFT the padded filter sequence of length "R" and save this complex 
> sequence.
> 
> 6) FFT each padded data sequence of length "R".
> 
> 7) Multiply the filter sequence by the data sequence (length "R").
> 
> 8) IFFT the resulting product (length "R").
> 
> The result will be a real time sequence which is nonzero for M + N -1 
> samples.
> 
> Note that N-1 nonzero samples at each end are filter transients.  This is 
> the same as if you applied the filter in time to a data sequence of length M 
> plus ending with enough zeros to see the tail end.
> 
> Fred 
> 
> 
thanks a lot for your introduction, depending on length of signal I 
reached factors between 40 and 110, as I guessed the old routine was not 
optimal for this problem.

Thomas

In article <gl1k28$ivp$1@ork-un.noris.net>, thomas.weber@logomotive.eu 
says...
>
>
>Hi,
>
>
>
>I've to apply some special shaped filter to long signals, the filter 
>coefficients I know already but the existing filter routine is probably to 
>slow. I want to update the filter routine without becoming a signal 
>processing expert ;-)
>
>
>
>What is the right way to apply a filter
>
>a)      in time domain or
>
>b)      in frequency domain
>
>
>
>It is possible to use some existing library code in FORTRAN ?

S/FILSYN will not only design just about any filter you can imagine but will 
automatically create a Fortran subroutine that implements the filter. (You 
call the subroutine once per sample.)

http://www.alkeng.com/sfilsyn.html

(Warning: It ain't cheap.)

In article <gl6mfd$8uq$1@ork-un.noris.net>, thomas.weber@logomotive.eu says...
>
>
>I've expected here some more useful information, like read this book, look 
>at this page or some example code.
>
>Things like
>
>    "Working in frequency domain is no big deal",
>
>    "Learning something new will take time" or
>
>    "there are a couple of key details you need to know to get things right"
>
>are not really helpful to start with the topic.
>
>
>
>I operate offline but if I apply one of the filter in our signal processing 
>software (NI Diadem) than it takes some seconds and in the existing Fortran 
>code some minutes (I have to apply this to many files and channels, so it is 
>hard to life with a factor in computation time of about 60), therefore I 
>look for a better filter code.
>
>BTW It don't have to be the best.

If the filter is FIR, then you can exploit the DOT_PRODUCT intrinsic in modern 
 Fortran (i.e., Fortran 90 or higher). In the time domain, an FIR filter is a 
convolution of the filter coefficients and the input data, so in Fortran, 
computing the filter's output sample is as simple as writing

output_sample=DOT_PRODUCT(coefficients(1:n),data(j:(j+n-1))

where coefficients is a real array containing n filter coefficients and data 
is a real array containing the last n input samples.

In Fortran, the cleanest way of updating the data array with each new sample 
is by using the EOSHIFT intrinsic, although it may be more computationally 
efficient to write a circular buffer explicitly.

The advantage of using the DOT_PRODUCT and EOSHIFT intrinsics is that they are 
probably highly optimized by the compiler.