MATLAB周波数領域でRam-Lakフィルタ（ランプフィルタ）を実装する方法は？

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は、カットオフ周波数であるように思わフィルターを使用しようと

は、私は、サンプリングレートとは何かを持っていますか？

また、2つの関数の畳み込みは、フーリエ空間における単純な乗算であると言われています。私は私があるとサイノグラムにこれを適用する方法がないアイデアを聞いていない、特に与えられていないBと、全くフィルタを実装する方法を理解していない

、私は誰かが私をここに助けることができると思います。私は2時間のグーグルで過ごし、ここで何が必要なのかを理解しようとしましたが、その実装方法を理解できませんでした。あなたが記載されている

Answer 1

式では、ラドンはフーリエ領域でのフィルタリングせずに逆変換をやってみたかった場合、中間結果です。別の方法としては、40年前より高速だったかもしれない空間領域での畳み込みを使用して、フィルタリングされた逆投影アルゴリズム全体を行うことです。最終的にあなたが投稿した数式を再現するでしょう。しかし、私は今はお勧めできません、特にあなたの最初の再建のためのものではありません。最初にヒルベルト変換を理解する必要があります。

とにかく、ここで必須のシェップ - ローガンファントムは、投影再構成をバックフィルタんいくつかのMatlabのコードです。私は、あなたがRam-Lakフィルタを使って独自のフィルタリングを行う方法を示します。私が本当に動機付けされていれば、ラドン/イラドンをいくつかのinterp2コマンドとサマリに置き換えます。

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