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

问题描述

我有一个任务来实现Ram-Lak过滤器，但是几乎没有给出任何信息（除了看fft，ifft，fftshift，ifftshift）。

我有一张我必须通过Ram-Lak过滤的正弦图。

我尝试使用过滤器

 <$ c如果I == 0 
 
（b ^ 2）/（2 * pi ^ 2）* 0，如果I even 
 
 -1 /（pi ^ 2 * I ^ 2）如果我奇数
  

b似乎是截止频率，我与采样率有关系吗？

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

我根本不知道如何实现滤波器，特别是没有给出，没有告诉我是什么，也不知道如何将其应用于正弦图，我希望有人能帮助我。我用了2个小时的搜索引擎，试图了解这里需要做什么，但是我不明白如何实现它。

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

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

$ b $ phantomData = phantom（）;

  N = size（phantomData，1）; 
 
 theta = 0：179; 
 N_theta =长度（theta）; 
 [R，xp] =氡（phantomData，theta）; 
 
％制作一个Ram-Lak过滤器。这只是abs（f）。 
 N1 =长度（xp）; 
 freqs = linspace（-1，1，N1）。 
 myFilter = abs（freqs）; 
 myFilter = repmat（myFilter，[1 N_theta]）; 
 
％做我自己的FT域过滤
 ft_R = fftshift（fft（R，[]，1），1）; 
 filteredProj = ft_R。* myFilter; 
 filteredProj = ifftshift（filteredProj，1）; 
 ift_R = real（ifft（filteredProj，[]，1））; 
 
％告诉matlab做没有过滤器的逆FBP 
 I1 = iradon（ift_R，theta，'linear'，'none'，1.0，N）; 
 
 subplot（1,3,1）; imagesc（real（I1））; title（'手动过滤'）
 colormap（灰色（256））;轴图像;轴偏离
 
％作为比较，请求matlab使用他们的Ram-Lak滤波器实现
 I2 = iradon（R，theta，'linear'，'Ram-Lak'，1.0，N） ; 
 
 subplot（1,3,2）; imagesc（real（I2））; title（'Matlab filtering'）
 colormap（gray（256））;轴图像;轴关
 
％为乐趣，重做过滤错误的目的
％排除高频创建一个低分辨率重建
 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（'低分辨率过滤'）
 colormap（灰色（256））;轴图像; axis off 
  

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.

解决方案

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

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