Posted by
Jacob Keller-2 on
Nov 01, 2018; 11:49pm
URL: http://imagej.273.s1.nabble.com/Frequency-Filtering-in-Time-Dimension-tp5021367p5021395.html
If you mainly want to suppress *one* frequency, you can use moving
> averages, with a kernel length equal to the period you want to subtract
> (i.e., a rectangular kernel; this gives you no oscillations). To do it
> in time (slice) direction, reslice the stack and use the 'Fast Filters'
> plugin. The kernel length is 2*radius+1. Maybe apply some weak Gaussian
> filter in addition.
>
My data are responses to a given stimulus frequency, so I want to extract
responses at just that frequency. I have actually been doing just what you
suggest with a plugin from Jay Unruh which subtracts a moving window
average--I set the window to the period, output both the moving window and
the subtracted stacks. It's not perfect, however, since the moving window
introduces its own artifacts, as one can see from Bode plots of moving
windows. I guess the ideal might be simply to fit the data with the right
equations/parameters, then plot those parameters over time. But that's not
as agnostic as a filter.
If you want a more square-like filter in the Fourier domain (and you
> don't care about ringing), and the desired cutoff frequency is rather
> high, you can use Process>Filter>Convolve on the resliced stack. For a
> perfectly sharp frequency cutoff you would need a sin(a*x)/(a*x) kernel
> (the so-called "sinc" function); but this would give you infinite kernel
> size.
>
https://en.wikipedia.org/wiki/Sinc_function> Multiply it with a Gaussian to limit the kernel size; this makes the
> cutoff a bit smoother (in frequency space, it gets convolved, i.e.,
> blurred, by the Fourier transform of the Gaussian that you have used for
> multiplying). The narrower the Gaussian you multiply with the sync, the
> smaller a kernel you can use (if you neglect very small numbers in the
> periphery of the kernel), but the frequency cutoff becomes smoother.
> Maybe a Gaussian with a sigma between the first and second zero is a
> good starting point (i.e., sigma between a*pi and a*2pi).
>
I was thinking of trying to use Lanczos filtering, which is two sinc's
superimposed. It's supposed to be a good compromise.
This is just an other way to express what Herbie has talked in his
> recent mails as "filtering in (real) time" vs. "filtering an image
> stack".
>
Okay, I get what Herbie was saying now. It would seem that the static case
would be much easier, then. I am thinking that fitting with some restraints
might be the best way to do this.
All the best,
Jacob
>
> Michael
> ________________________________________________________________
>
>
> On 2018-10-31 17:32, Jacob Keller wrote:
> > Thanks Michael--this is really clear and helpful! I checked out some
> > Bode
> > plots of various filters, and it seems Gaussian is pretty non-brick
> > wall. I
> > wonder whether a Butterworth or Bessel filter has been implemented?
> >
> > JPK
> >
> > On Wed, Oct 31, 2018 at 7:36 AM Michael Schmid
> > <
[hidden email]>
> > wrote:
> >
> >> Hi Jacob,
> >>
> >> when you are doing a Gaussian Blur, in the Fourier domain it
> >> corresponds
> >> to multiplying the amplitudes with a Gaussian as well.
> >>
> >> A Gaussian
> >> exp(-pi * x²/a²)
> >> in real space corresponds to
> >> exp(-pi * f² * a²)
> >> in Fourier space
> >>
> >> If you want to attenuate a frequency component f1 by a factor 1/e =
> >> 0.37, you thus need a² = 1/(pi f1²), i.e. exp (pi² * f1² * x²)
> >>
> >> In terms of sigma, a Gaussian is given by
> >> exp(-x²/(2*sigma²))
> >>
> >> Thus, attenuation by a factor of 1/e at f1 corresponds to
> >> sigma = 1/(pi * f1 * sqrt(2))
> >>
> >>
> >> You can easily check it with the following macro:
> >>
> >> period = 16; // between about 10 and 100
> >> // best accuracy with powers of 2
> >> f=1/period; // spatial frequency
> >> newImage("Untitled", "32-bit white", 256, 256, 1);
> >> run("Macro...", "code=v=sin(2*PI*x*"+f+")"); //create sine wave
> >> sigma = 1/(PI*f*sqrt(2));
> >> run("Gaussian Blur...", "sigma="+sigma);
> >> run("Set Measurements...", "min display redirect=None decimal=4");
> >> makeRectangle(25, 0, 200, 256);
> >> run("Measure");
> >>
> >>
> >> You will see that the maxima and minima of the sine wave (which has
> >> had
> >> an original amplitude of 1) gets attenuated to an amplitude of about
> >> 0.37.
> >>
> >> If you take sigma = 1/(pi * f1) you will get an attenuation factor of
> >> 1/e² = 0.13 at the frequency f1.
> >>
> >> So far an excursion into the math of Gaussians and their Fourier
> >> transform...
> >>
> >>
> >> Michael
> >> ________________________________________________________________
> >> On 30.10.18 16:23, Jacob Keller wrote:
> >> >>
> >> >> do you have simply a 1D data set or an image stack, with slices
> >> for
> >> >> different times?
> >> >>
> >> >
> >> > The latter--image stack, one plane over time.
> >> >
> >> > - In the second case, you can use Process>Filters>Gaussian Blur 3D
> >> and
> >> >> specify the x and y sigma as zero.
> >> >> For high-pass filtering, you would have to duplicate the stack
> >> first
> >> and
> >> >> then subtract the filtered image.
> >> >>
> >> >
> >> > Ah, this is a great idea--I will try it out. How can I figure out
> >> the
> >> > relationship between sigma and the desired frequency cutoff? For
> >> example,
> >> > let's say the stack has 120 frames per period--what sigma value
> >> should
> >> be
> >> > input?
> >> >
> >> > Thanks very much for these suggestions--I think they are going to
> >> be
> >> > very helpful,
> >> >
> >> > Jacob
> >> >
> >> > If you rather want moving averages or a median, you can reslice the
> >> >> stack (Image>Stacks>Reslice) to have the time direction in x or y
> >> and
> >> >> apply the 'Fast Filters' plugin, which can do 1D filtering. Then
> >> Reslice
> >> >> to make time the z axis again.
> >> >> The 'Fast Filters' plugin also has an option to subtract the
> >> filtered
> >> >> image (i.e., highpass operation)
> >> >>
> >> >>
> >> >> Michael
> >> >> ________________________________________________________________
> >> >> On 29.10.18 20:51, Jacob Keller wrote:
> >> >>> Hi All,
> >> >>>
> >> >>> is there a way to do high-pass/low-pass filtering in the time
> >> dimension?
> >> >> Or
> >> >>> a suggestion for a workaround?
> >> >>>
> >> >>> All the best,
> >> >>>
> >> >>> Jacob Keller
> >> >>
> >> >> --
> >> >> ImageJ mailing list:
http://imagej.nih.gov/ij/list.html> >> >>
> >> >
> >> > --
> >> > ImageJ mailing list:
http://imagej.nih.gov/ij/list.html> >> >
> >>
> >> --
> >> ImageJ mailing list:
http://imagej.nih.gov/ij/list.html> >>
> >
> > --
> > ImageJ mailing list:
http://imagej.nih.gov/ij/list.html>
> --
> ImageJ mailing list:
http://imagej.nih.gov/ij/list.html>
--
ImageJ mailing list:
http://imagej.nih.gov/ij/list.html