MATLAB如何在频域实现Ram-Lak滤波器（Ramp滤波器）？

发布于 2024-11-24 08:13:05 字数 482 浏览 0 评论 0原文

我有一个实现 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好像是截止频率，我和采样率有关系吗？

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

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

原文

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

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

I try to use the filter

                      1/4              if I == 0

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

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

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.

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.

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

圈圈圆圆圈圈 2024-12-01 08:13:05

如果您想进行 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

手动过滤演示

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.

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();

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

Demonstration of manual filtering

回复收藏 0 原文

~没有更多了~