MATLAB 如何在频域中实现 Ram-Lak 滤波器(Ramp filter)? [英] MATLAB How to implement a Ram-Lak filter (Ramp filter) in the frequency domain?

问题描述

我有一个任务来实现一个 Ram-Lak 过滤器，但几乎没有给出关于它的信息(除了看 fft、ifft、fftshift、ifftshift).

I have an assignment to implement a Ram-Lak filter, but nearly no information given on it (except look at fft, ifft, fftshift, ifftshift).

我有一个正弦图，我必须通过 Ram-Lak 进行过滤.还给出了投影的数量.

I have a sinogram that I have to filter via Ram-Lak. Also the number of projections is given.

我尝试使用过滤器

                      1/4              if I == 0

(b^2)/(2*pi^2)  *     0                if I even

                      -1/(pi^2 * I^2)  if I odd

b好像是截止频率，跟我采样率有关系吗?

b seems to be the cut-off frequency, I has something to do with the sampling rate?

也有人说两个函数的卷积是傅立叶空间中的简单乘法.

Also it is said that the convolution of two functions is a simple multiplication in Fourier space.

我根本不明白如何实现过滤器，尤其是没有给出 b，没有告诉我是什么，也不知道如何将它应用于正弦图，我希望有人能在这里帮助我.我花了 2 小时的谷歌搜索并试图了解这里需要做什么，但我无法理解如何实现它.

I do not understand how to implement the filter at all, especially with no b given, not told what I is and no idea how to apply this to the sinogram, I hope someone can help me here. I spent 2hrs googling and trying to understand what is needed to do here, but I could not understand how to implement it.

推荐答案

如果您想在傅立叶域中进行不带滤波的逆 Radon 变换，则您列出的公式是一个中间结果.另一种方法是在空间域中使用卷积来完成整个滤波反投影算法，这在 40 年前可能会更快；你最终会重新推导你发布的公式.但是，我现在不推荐它，尤其是对于您的第一次重建；你应该先真正了解希尔伯特变换.

The formula you listed is an intermediate result if you wanted to do an inverse Radon transform without filtering in the Fourier domain. An alternative is to do the entire filtered back projection algorithm using convolution in the spatial domain, which might have been faster 40 years ago; you would eventually rederive the formula you posted. However, I wouldn't recommended it now, especially not for your first reconstruction; you should really understand the Hilbert transform first.

无论如何，这里有一些 Matlab 代码，它执行强制性的 Shepp-Logan 幻像滤波反投影重建.我将展示如何使用 Ram-Lak 过滤器进行自己的过滤.如果我真的有动力，我会用一些 interp2 命令和求和来替换 radon/iradon.

Anyway, here's some Matlab code which does the obligatory Shepp-Logan phantom filtered back projection reconstruction. I show how you can do your own filtering with the Ram-Lak filter. If I was really motivated, I would replace radon/iradon with some interp2 commands and summations.

<代码>phantomData=phantom();

phantomData=phantom();

N=size(phantomData,1);

theta = 0:179;
N_theta = length(theta);
[R,xp] = radon(phantomData,theta);

% make a Ram-Lak filter. it's just abs(f).
N1 = length(xp);
freqs=linspace(-1, 1, N1).';
myFilter = abs( freqs );
myFilter = repmat(myFilter, [1 N_theta]);

% do my own FT domain filtering
ft_R = fftshift(fft(R,[],1),1);
filteredProj = ft_R .* myFilter;
filteredProj = ifftshift(filteredProj,1);
ift_R = real(ifft(filteredProj,[],1));

% tell matlab to do inverse FBP without a filter
I1 = iradon(ift_R, theta, 'linear', 'none', 1.0, N);

subplot(1,3,1);imagesc( real(I1) ); title('Manual filtering')
colormap(gray(256)); axis image; axis off

% for comparison, ask matlab to use their Ram-Lak filter implementation
I2 = iradon(R, theta, 'linear', 'Ram-Lak', 1.0, N);

subplot(1,3,2);imagesc( real(I2) ); title('Matlab filtering')
colormap(gray(256)); axis image; axis off

% for fun, redo the filtering wrong on purpose
% exclude high frequencies to create a low-resolution reconstruction
myFilter( myFilter > 0.1 ) = 0;
ift_R = real(ifft(ifftshift(ft_R .* myFilter,1),[],1));
I3 = iradon(ift_R, theta, 'linear', 'none', 1.0, N);

subplot(1,3,3);imagesc( real(I3) ); title('Low resolution filtering')
colormap(gray(256)); axis image; axis off

这篇关于MATLAB 如何在频域中实现 Ram-Lak 滤波器(Ramp filter)?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！