1D ramp filter
Alistair Templeton
I'm playing with iradon, but would like to try filtering the sinogram first but I can't figure out how to design the 1D filter to apply to the columns.
If I just have a 1-D array of numbers, how would I apply a simple ramp filter to it?
Matt J
> Way simple question, it's totally stumped me.
>
> I'm playing with iradon, but would like to try filtering the sinogram first but I can't figure out how to design the 1D filter to apply to the columns.
>
> If I just have a 1-D array of numbers, how would I apply a simple ramp filter to it?
help conv
help fft
help filter
Alistair Templeton
> help conv
> help fft
> help filter
apologies, I should have indicated that I've found my way to the above commands but am having trouble getting through.
I've tried the following:
a = [ 1 5 6 2 6 3 4 5]; %nonsense sample data
fa = fft(a);
ramp = (1:8) / 8; % what i'm hoping is a basic ramp filter.
filtered_a = ifft(fa.*ramp); % gives me the wrong result
This isn't homework, it's just a side thing I'd like to add to a presentation so I'd very much appreciated being given the fish as I know this is trivial for someone who knows how to do it!
cheers
Matt J
> apologies, I should have indicated that I've found my way to the above commands but am having trouble getting through.
>
> I've tried the following:
>
> a = [ 1 5 6 2 6 3 4 5]; %nonsense sample data
> fa = fft(a);
> ramp = (1:8) / 8; % what i'm hoping is a basic ramp filter.
> filtered_a = ifft(fa.*ramp); % gives me the wrong result
===============
This is only a 1-sided ramp. What you really would need is something like
ramp=abs((1:8)-5);
Incidentally, even though both you in the above and also iradon design the filter kernel in the frequency domain, it is a bad idea to do so.
There are closed form formulas for the ramp and other common filter kernels in the non-Fourier domain. It is therefore recommendable to build the filter kernel in the non-Fourier domain first and then take its FFT later for the purpose of executing the convolution efficiently. This greatly cuts down on aliasing and the need for zero-padding.
For details, see Section 3.3.3 in
@BOOK{
kak:88,
author = {A. C. Kak and M. Slaney},
title = {Principles of computerized tomographic imaging},
publisher = {IEEE Press},
address = {New York},
year = 1988
}
Matt J
> This is only a 1-sided ramp. What you really would need is something like
>
> ramp=abs((1:8)-5);
Actually, you would want
ramp=ifftshift( abs((1:8)-5) );
so as to deal with FFT circulancy issues.
Alistair Templeton
>
> ramp=ifftshift( abs((1:8)-5) );
>
> so as to deal with FFT circulancy issues.
that's exactly what I was missing, thank you!
I was googling around in the imaging crowd and apparently everybody and their mother, in that circle, designs their filters in the frequency domain. I suppose it's much easier to visualize that way.
Matt J
> I was googling around in the imaging crowd and apparently everybody and their mother, in that circle, designs their filters in the frequency domain. I suppose it's much easier to visualize that way.
==================
It's an easy and common mistake to make, but that doesn't mean everyone is doing it.
Maybe everyone's mothers is, but not necessarily their sisters.