71516

Адаптивная фильтрация двумерных сигналов

Лабораторная работа

Коммуникация, связь, радиоэлектроника и цифровые приборы

В отличие от imnoise, М-функция imnoise2 порождает шумовую матрицу R размера MxN., которая не нормируется. Другое значительное отличие от функции imnoise состоит в том, что выходом imnoise служит зашумленное изображение, a imnoise2 порождает только шумовую матрицу.

Русский

2014-11-08

3.57 MB

8 чел.

МИНОБРНАУКИ РОССИИ

Государственное Образовательное учреждение

Высшего Профессионального Образования

ПЕНЗЕНСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНОЛОГИЧЕСКИЙ УНИВЕРСИТЕТ

Кафедра «Информационные технологии и системы»

Дисциплина «Цифровая обработка сигналов»

Отчет по лабораторной работе №3

«Адаптивная фильтрация двумерных сигналов»

Выполнила: студентка группы 09И

Соляникова Е.А.

Проверил: к.т.н., доцент кафедры

Кирин Ю.П.

Пенза, 2013 г.

Цель работы:

  1.  Освоение принципов адаптивной фильтрации в частотной области.
  2.  Овладение навыками адаптивной фильтрации в частотной области.

Выполнение:

В отличие от imnoise, М-функция imnoise2 порождает шумовую матрицу R размера MxN., которая не нормируется. Другое значительное отличие от функции imnoise состоит в том, что выходом imnoise служит зашумленное изображение, a imnoise2 порождает только шумовую матрицу. Пользователь должен сам задавать желаемые параметры шума. Заметьте, что шумовая матрица, генерируемая по методу «соль и перец», принимает три возможных значения: 0 - это «перец», 1 - «соль», а 0.5 соответствует отсутствию шума. Эту матрицу необходимо еще обрабатывать для дальнейшего использования. Например, чтобы испортить изображение этим шумом, необходимо сначала найти (с помощью функции find) все координаты шума R, равные 0, и присвоить соответствующим пикселям изображения самое малое значение (обычно это 0). Затем следует найти все координаты шума R, равные 1, и присвоить соответствующим пикселям изображения самое большое допустимое значение (обычно это 255 для 8-ми битовых изображений). Эта процедура моделирует практическое воздействие шума «соль и перец» на изображение. Ниже приведен листинг функции imnoise2: 

function R=imnoise2(type, M, N, a, b)

if nargin ==1

   a=0;

   b=1;

   M=1; N=1;

elseif nargin ==3

   a=0;

   b=1;

end

switch lower(type)

   case 'uniform'

       R=a+(b-a)*rand(M,N);

   case 'gaussian'

       R=a+b*randn(M,N);

   case 'salt & pepper'

       if nargin <=3

           a=0.05; b=0.05;

       end

       

       if (a+b)>1

           error('The sum Pa + Pb is not > 1')

       end

       R(1:M, 1:N)=0.5;

       X=rand(M,N);

       c=find(X<=a);

       R(c)=0;

       u=a+b;

       c=find(X>a & X<=u);

       R(c)=1;

   case 'lognormal'

       if nargin <=3

           a=1; b=0.25;

       end

       R=exp(b*randn(M,N)+a);

   case 'rayleigh'

       R=a+(-b*log(1-rand(M,N))).^0.5;

   case 'exponential'

       if nargin <= 3

           a=1;

       end

       if a <= 0

           error('Parameter a must be positive for exponential type')

       end

       k=-1/a;

       R=k*log(1-rand(M,N));

   case 'erlang'

       if nargin <=3

           a=2; b=5;

       end

       if (b~=round(b) | b<=0)

          error('Param b must be a positive integer for Erlang')

       end

       k=-1/a;

       R=zeros(M,N);

       for j=1:b

           R=R+k*log(1-rand(M,N));

       end

   otherwise

       error('Unknown distribution type')

end

Следующая М-функция допускает произвольное число местоположений импульсов (координат частот), каждый со своей амплитудой, частотой и сдвигом фазы, и она вычисляет r(х, у) в виде суммы синусоид в форме:

function [r, R, S]=imnoise3(M, N, C, A, B)

[K,n]= size(C);

if nargin ==3

   A(1:K)=1.0;

   B(1:K, 1:2)=0;

elseif nargin ==4

   B(1:K, 1:2)=0;

end

R=zeros(M,N);

for j=1:K

   u1=M/2+1+C(j,1);

   v1=N/2+1+C(j,2);

   R(u1,v1)=i*(A(j)/2)*exp(i*2*pi*C(j,1)*B(j,1)/M);

   u2=M/2+1-C(j,1);

   v2=N/2+1+C(j,2);

   R(u2,v2)=-i*(A(j)/2)*exp(i*2*pi*C(j,2)*B(j,2)/N);

end

S=abs(R);

r=real(ifft2(ifftshift(R)));

end

На рис. 1 и рис. 2 построен спектр и пространственный синусовый шум с помощью последовательности команд:

C=[0 64; 0 128; 32 32; 64 0; 128 0; -32 32];

[r, R, S]=imnoise3(512, 512, C);

figure, imshow(S, [])

figure, imshow(r, [])

На рис. 3 и рис. 4 показан результат повторения предыдущих команд, но с параметрами:

C=[0 32; 0 64; 16 16; 32 0; 64 0; -16 16];

[r, R, S]=imnoise3(512, 512, C);

figure, imshow(S, [])

figure, imshow(r, [])

Аналогично, рисунок 5 построен при:

C=[6 32; -2 2];

[r, R, S]=imnoise3(512, 512, C);

figure, imshow(r, [])

Рисунок 6 получен с тем же вектором С, но с модифицированным вектором амплитуд:

C=[6 32; -2 2];

A=[1 5];

[r, R, S]=imnoise3(512, 512, C, A);

figure, imshow(r, [])

Из рисунка 7 видно, что на изображении доминируют синусовые волны низкой частоты. Это легко понять, так как их амплитуда в пять раз больше амплитуды высокочастотной компоненты.

Результат выполнения программ:

 

Рис. 1 - Спектр заданных импульсов Рис. 2 - Соответствующий синусовый шум

 

Рис. 3 - Спектр заданных импульсов Рис. 4 - Соответствующий синусовый шум

 

Рис. 5 – Другой шумовой образ Рис. 6 – Другой шумовой образ

  1.  Использование функции spfilt, которая совершает фильтрацию в пространственной области с помощью любого фильтра.

Ниже приведен текст функции spfilt:

function f=spfilt(g, type, m, n, parameter)

if nargin ==2

   n=3; m=3; q=1.5; d=2;

elseif nargin==5

   q=parameter; d=parameter;

elseif nargin ==4

   q=1.5; d=2;

else

   error('Неверные входные параметры')

end

switch type

   case 'amean'

       w=fspecial('average', [m,n]);

       f=imfilter(g, w, 'replicate');

   case 'gmean'

       f=gmean(g,m,n);

   case 'hmean'

       f=harmean(g,m,n);

   case 'chmean'

       f=charmean(g,m,n,q);

   case 'median'

       g=medfilt2(g, [m,n], 'symmetric');

   case 'max'

       f=ordfilt2(g, m*n, ones(m,n)*0.4, 'symmetric');

   case 'min'

       f=ordfilt2(g, 1, ones(m,n)*0.4, 'symmetric');  

   case 'midpoint'

       f1=ordfilt2(g, m*n, ones(m,n)*0.4, 'symmetric');

       f2=ordfilt2(g, m*n, ones(m,n)*0.4, 'symmetric');

       f=imlincomb(0.5, f1, 0.5, f2);

   case 'atrimmed'

       if (d<0)|(d/2~=round(d/2))

           error('Параметр d должен быть положительным целым')

       end

       f=alphatrim(g,m,n,d);

   otherwise

       error('Неизвестный тип фильтра')

end

___________________________________________________________

function f=gmean(g,m,n)

inclass=class(g);

g=im2double(g);

warning off;

f=exp(imfilter(log(g), ones(m,n), 'replicate')).^(1/m/n);

warning on;

end

___________________________________________________________

function f=harmean(g,m,n,q)

inclass=class(g);

g=im2double(g);

f=m*n./imfilter(1./(g+eps), ones(m,n), 'replicate');

end

___________________________________________________________

function f=charmean(g,m,n,q)

inclass=class(g);

g=im2double(g);

f=imfilter(g.^(q+1), ones(m,n), 'replicate');

f=f./(imfilter(g.^q, ones(m,n), 'replicate')+eps);

end

___________________________________________________________

function f=alphatrim(g,m,n,d)

inclass=class(g);

g=im2double(g);

f=imfilter(g, ones(m,n), 'symmetric');

for k=1:d/2

f=imsubstract(f, ordfilt2(g, k, ones(m,n), 'symmetric'));

end

for k=(m*n-(d/2)+1):m*n

f=imsubstract(f, ordfilt2(g, k, ones(m,n), 'symmetric'));

end

f=f/(m*n-d);

end

На рисунке 7 приведено изображение класса uint8, испорченное только шумом «перец» с вероятностью 0.2. Это изображение построено следующими командами:

f=imread('beaut.jpg');

[M,N]=size(f);

R=imnoise2('salt & pepper', M, N, 0.2, 0);

c = find(R == 0);

gp = f;

gp(c) = 0;

figure, imshow(gp)

Рисунок 8, подпорченный только шумом «соль», получен командами:

R = imnoise2('salt & pepper', M, N, 0, 0.2);

c = find(R == 1);

gs = f;

gs(c) = 255;

figure, imshow(gp)

Избавиться от шума «перец» можно с помощью контргармонического фильтра с положительным значением Q. Рисунок 9 построен командами:

fp = spfilt(gp, 'chmean', 3, 3, 1.5);

figure, imshow(fp)

Аналогично для борьбы с шумом «соль» можно воспользоваться контргармоническим фильтром с отрицательным значением Q:

fs = spfilt(gs, 'chmean', 3, 3, -1.5);

figure, imshow(fs)

На рисунке 10 показан результат. Можно также действовать минимальным и максимальным фильтрами. Например, рисунки 11 и рис. 12 получены, соответственно, из рис. 7 и рис. 8 с помощью следующих команд:

fpmax = spfilt(rgb2gray(gp), 'max', 3, 3);

figure, imshow(fpmax)

fpmin = spfilt(rgb2gray(gs), 'min', 3, 3);

figure, imshow(fpmin)

Результат выполнения программ:

 

Рис. 7 - Изображение, испорченное  Рис. 8 - Изображение, испорченное

шумом «перец» с вероятностью 0.2 шумом «соль» с вероятностью 0.2

 

Рис. 9 - Результат фильтрации Рис. 10 - Контргармоническим фильтром размера 3х3 порядка Q=1.5

 

Рис. 11 - Результат фильтрации Рис. 12 - Контргармоническим фильтром при Q=-1.5

  1.  Адаптивная медианная фильтрация:

В качестве примера рассмотрим изображение f на рис. 7, испорченный шумом «соль и перец», который наложен на изображение командой:

g=imnoise(f, 'salt & pepper', 0.25);

figure, imshow(g)

а на рисунке 13 приведен результат выполнения следующей команды:

f1 = medfilt2(rgb2gray(g), [7 7], 'symmetric');

figure, imshow(f1)

Это изображение, очевидно, очищено от шума, но одновременно с этим оно стало слишком расплывчатым и искаженным. Напротив, совершив действие:

f2 = adpmedian(rgb2gray(g), 7);

figure, imshow(f2)

мы получим изображение 14, на котором шум также отсутствует, но оно существенно четче, чем рис. 15.

Результат выполнения программы:

 

Рис. 13 - Изображение, испорченное  Рис. 14 - Результат применения

шумом «соль и перец» с плотностью  медианного фильтра размером 7х7

ошибок 0.25

Рис. 15 - Результат адаптивной медианной фильтрации с Smax=7

Ниже приведен текст функции adpmedian.

function f=adpmedian(g, Smax)

if (Smax<=1)|(Smax/2==round(Smax/2))|(Smax~=round(Smax))

   error('SMAX must be an odd integer > 1')

end

[M,N]=size(g);

 

f=g;

f(:)=0;

alreadyProcessed=false(size(g));

 

for k=3:2:Smax

   zmin=ordfilt2(g,1,ones(k,k), 'symmetric');

   zmax=ordfilt2(g,k*k,ones(k,k), 'symmetric');

   zmed=medfilt2(g,[k k], 'symmetric');

   

   processUsingLevelB=(zmed > zmin)&(zmax > zmed) & ...

       ~alreadyProcessed;

   zB=(g > zmin)&(zmax > g);

   outputZxy = processUsingLevelB & zB;

   outputZmed = processUsingLevelB & ~zB;

   f(outputZxy)=g(outputZxy);

   f(outputZmed)=zmed(outputZmed);

   

  alreadyProcessed=alreadyProcessed | processUsingLevelB;

  if all(alreadyProcessed(:))

      break;

  end

end

f(~alreadyProcessed)=zmed(~alreadyProcessed);

end

  1.  Моделирование размытого зашумленного изображения:

Многие алгоритмы восстановления изображений работают медленно с большими изображениями. Поэтому имеет смысл экспериментировать с малыми изображениями для сокращения времени вычислений, что позволяет скорее оптимизировать алгоритм. В этом случае для целей визуализации полезно иметь возможность увеличивать изображение посредством дублирования пикселов. Это делается с помощью следующей функции:

function B=pixeldup(A,m,n)

if nargin < 2

   error('At least two inputs are required')

end

if nargin ==2

   n=m;

end

u=1:size(A,1);

m=round(m);

u=u(ones(1,m),:);

u=u(:);

 

v=1:size(A,2);

n=round(n);

v=v(ones(1,n),:);

v=v(:);

B=A(u,v);

End

На рисунке 16 построено изображение с помощью команды:

f=checkerboard(8);

figure, imshow(f)

Испорченное изображение 17 получено командами:

PSF = fspecial('motion', 7, 45);

gb=imfilter(f, PSF, 'circular');

figure, imshow(gb)

Шумовое изображение на рисунке 18 сгенерировано командами:

noise=imnoise(zeros(size(f)), 'salt & pepper', 0.2);

gb=gb+noise;

figure, imshow(gb)

Результат выполнения программы:

 

Рис. 16 – Исходное изображение  Рис. 17 - Изображение, размытое функцией fspecial с параметрами len=7 и  theta = 45

Рис. 18 - Зашумленное изображение

  1.  Винеровская фильтрация:

Применение функции deconvwnr при восстановлении размытого зашумленного изображения. Рисунок 19 совпадает с рис. 18, а рис. 20 получен командой:

fr1=deconvwnr(gb, PSF);

figure, imshow(fr1)

где g - это искаженное изображение, a PSF - функция разброса точек. Как уже отмечалось выше, frl - это результат прямого применения инверсной фильтрации, и, как ожидалось, на этом изображении преобладают шумовые эффекты.

Число R было получено с помощью исходного и шумового изображений по формулам:

Sn=abs(fft2(noise)).^2;

nA=sum(Sn(:))/prod(size(noise));

Sf=abs(fft2(f)).^2;

fA=sum(Sf(:))/prod(size(noise));

R=nA/fA;

Для восстановления изображения с помощью этого соотношения R выполним команду:

fr2=deconvwnr(gb, PSF, R);

figure, imshow(fr2)

Из рисунка 21 видно, что такой подход дает существенное улучшение по сравнению с инверсной фильтрацией. Наконец, используем автокорреляционные функции при восстановлении изображения (напомним о необходимости использования функции fftshift при центровании):

NCDRR=fftshift(real(ifft2(Sn)));

ICDRR=fftshift(real(ifft2(Sf)));

fr3=deconvwnr(gb, PSF, NCDRR, ICDRR);

figure, imshow(fr3)

Из рисунка 22 видно, что результат этих операций весьма близок к оригиналу, хотя небольшой шум все еще проявляется.

Результат выполнения программы:

 

Рис. 19 - Размытое и зашумленное   Рис. 20 - Результат инверсной

изображение     фильтрации

 

Рис. 21 - Результат винеровской   Рис. 22 - Результат винеровской

фильтрации с постоянным    фильтрации с использованием

соотношением сигнал/шум   автокорреляционных функций

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


 

А также другие работы, которые могут Вас заинтересовать

52699. Повторение изученного материала. (Сложение и вычитание многозначных чисел) 51 KB
  Проверка выполнения домашнего задания. В течение урока мы будем помогать строителям дома выполняя разные математические задания тем самым будем закреплять изученный материал. а Подготовка площадки для строительства нового дома: 1.10кг б как вы думаете с чего начинается строительство дома ответы детей с планапроекта грамотное составление плана проекта способствует правильному строительству дома поэтому внимание: расшифруйте фразу из ответов решив эти примеры по 2 командам.