Butterworth filter using real coefficients
Rick Danielson
06-10-2008, 11:27 AM
Hi All,
Thanks to NR, I have the real and imaginary Fourier components of some real-valued data (sea surface temperature) and my fifth-order lowpass "Butterworth" filter is
L(f) = [1 + (f / f_c)^10]^0.5
where f is frequency and f_c is the cutoff. Of course, this isn't your everyday Butterworth filter because it's real-valued. (I apply it to both real and imaginary parts of my spectrum.) If I use the more common imaginary-valued filter, that will lead to phase changes in my filtered output, and I don't want that. One solution is to filter twice (once forward and once backward - the filtfilt function in matlab). However, it seems simpler just to use this real-valued filter above. Does anyone see the problem here? I'm posting this because I can't find any discussion about the use of real-valued filters like the one above - either positive or negative comments - besides what is given in NR. (Comments/citations are both welcome in this thread!)
Rick Danielson
07-10-2008, 11:57 AM
Perhaps the question wasn't as well posed as it could have been (but feedback is still welcome). I should have written L(f) = [1 + (f / f_c)^10]^-1 and called this the modulus squared of a fifth-order Butterworth filter. It is real-valued (multiplies both the real and imaginary spectral coefficients), with the desired result that the filter is zero-phase. The closest external reference to it that I've found is by Gustafsson (1996, Determining the initial states in forward-backward filtering, Trans. Signal Process., 46, 988--992). Interestingly, this is one of the references given for forward-backward filtering in Matlab, but there's no direct comparison between spectral and temporal filtering (although the advice in NR to stick to spectral sounds good to me)...
davekw7x
07-10-2008, 09:01 PM
Perhaps the question...
I think I understood the question OK; I just didn't know (and still don't know) what you would like to hear from us.
You can't prove something works by testing it, but sometimes you can get some insight (and sometimes you can determine that certain things don't work the way that you thought they did) by actually trying things.
It seems to me that, if you have a project involving digital filtering, you would investigate the different approaches by trying them. I mean, try them with some special test data for which you can know the "right" answer and is typical of what you will be using.
So: Have you tried any of the above? What are the problems with each approach? Can a zero-phase non-constant function of frequency be causal? What are the effects of approximating the desired response by a causal process in each case? Under what circumstances can these effects be taken into account in such a way as to make the approximation acceptable in the context of actual requirements?
Since you are considering the "read it forward and again backward" are we to think that you have the entire data set at your disposal (rather than wanting to filter it on the fly and getting an running output in the time domain)? And that you are going to perform a DFT on the entire data set (and then do the inverse DFT)? Are computational requirements (CPU time, real time, etc.) a factor in your selection of a method?
It seems to me that the approach using the Matlab filtfilt function that you mention would be very straightforward. You could also use the Matlab fft function, multiply by your filter equation and do the inverse with ifft so that you can compare the methods. So, what's the problem?
Of course there are other approaches.
What are your actual requirements? (Filter cutoff frequency relative to sampling frequency; maximum attenuation in the pass band; minimum attenuation in the cutoff band; stuff like that)? I mean, in principle, constant group delays can be attained with FIR filters that can be constructed to whatever response requirements that you specify. The delay (time from input samples to output results) will be about half the length of the filter, and there are no limits to the length of the data.
So, what kind of help, exactly, are you requesting?
Regards,
Dave