如何在频域中实现Ram-Lak滤波器（斜坡滤波器）的MATLAB代码？

Question

我有一个任务要实现Ram-Lak滤波器，但几乎没有提供任何信息（除了看fft、ifft、fftshift、ifftshift）。

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

我尝试使用这个滤波器。

                      1/4              if I == 0

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

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

“b”似乎是截止频率，而“I”与采样率有关？

此外，据说两个函数的卷积在傅里叶空间中是简单的乘法。

我完全不理解如何实现滤波器，特别是没有给出“b”，没有告诉我“I”是什么，也不知道如何将其应用于正弦图像，希望有人能够在这里帮助我。我花了2个小时搜索和尝试理解需要做些什么，但我无法理解如何实现它。

- mmlac

---

你的问题非常不清楚。I是什么？我们如何找到b？我查看的所有文本都没有使用相同的符号。Ram Lak滤波器只是一个abs(f)函数，就像双斜坡一样。如果你告诉我这些变量是什么，我就能帮助你了。 - Phonon

我的问题是我不知道。在文献中（即“非衍射源重建算法”第72页-链接），他们只使用了一个I（或teta，在那里是采样间隔）。 您能帮我基于正弦图和投影数量实现一个简单的abs()过滤器吗？ - mmlac

试试这个链接。http://laskin.mis.hiroshima-u.ac.jp/Kougi/07a/PIP/PIP12pr.pdf 这是一篇很好的关于它的文章。图2B是一个Ram Lak滤波器。 - Phonon

Answer 1

你列出的公式是在傅里叶域中不进行过滤进行反Radon变换的中间结果。另一种选择是使用空间域中的卷积执行整个滤波反投影算法，这在40年前可能会更快；你最终将重新推导出你发帖的公式。无论如何，这里有一些Matlab代码，可以进行基本的Shepp-Logan幻像滤波反投影重建。我展示了如何使用Ram-Lak滤波器进行自定义滤波。如果我真的有动力，我会用一些interp2命令和求和取代radon / iradon。 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

- user244795