MATLAB Как реализовать фильтр Рам-Лака (Ramp-фильтр) в частотной области? - PullRequest
Новая
       31

MATLAB Как реализовать фильтр Рам-Лака (Ramp-фильтр) в частотной области?

4 голосов
/ 15 июля 2011

У меня есть задание для реализации фильтра Рам-Лака, но информация о нем почти не указана (за исключением просмотра fft, ifft, fftshift, ifftshift).

У меня есть синограмма, которую я должен отфильтровать через Рам-Лака. Также приводится количество прогнозов.

Я пытаюсь использовать фильтр 

                      1/4              if I == 0

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

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

b кажется частотой среза, я как-то связан с частотой дискретизации?

Также сказано, что свертка двух функций является простым умножением в пространстве Фурье.

Я вообще не понимаю, как реализовать фильтр, особенно без б, не сказал, кто я, и не знаю, как применить это к синограмме, я надеюсь, что кто-то может помочь мне здесь. Я потратил 2 часа на поиски в Интернете и пытаясь понять, что нужно делать здесь, но я не мог понять, как это реализовать.

1 Ответ

10 голосов
/ 23 июля 2011

Формула, которую вы перечислили, является промежуточным результатом, если вы хотите выполнить обратное преобразование Радона без фильтрации в области Фурье.Альтернативой является создание всего алгоритма обратной проекции с фильтрацией с использованием свертки в пространственной области, что могло бы быть быстрее 40 лет назад;в конечном итоге вы получите новую формулу.Однако я бы не рекомендовал это сейчас, особенно не для вашей первой реконструкции;сначала вы должны по-настоящему понять преобразование Гильберта.

В любом случае, вот некоторый код Matlab, который выполняет обязательную реконструкцию обратной проекции с фильтром Фантома Шеппа-Логана.Я покажу, как вы можете сделать свою собственную фильтрацию с помощью фильтра Рам-Лака.Если бы я был действительно мотивирован, я бы заменил радон / ирадон некоторыми командами и суммами interp2.

phantomData=phantom();</p> <pre><code>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

...