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 - Результат винеровской

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

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

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


 

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

32875. АНТИЧНАЯ ФИЛОСОФИЯ. ПИФАГОР, СОКРАТ, ПЛАТОН, АРИСТОТЕЛЬ, ДЕМОКРИТ, ЭПИКУР 64.08 KB
  Подобно тому как математик ставит вопрос что такое единица и дает довольнотаки сложное определение этого казалось бы простейшего понятия так и философ с глубокой древности задается проблемой: что такое бытие что значит быть Эта специфика философии проливает известный свет и на вопрос о том почему и когда философия возникает. Философ ставит все это под сомнение но делает это для того чтобы докопаться до подлинных корней из которых растет сама данная традиция; в этом и состоит смысл его вопроса: что значит быть что такое бытие А...
32876. СРЕДНЕВЕКОВАЯ ФИЛОСОФИЯ. ОТЦЫ ЦЕРКВИ. ВАСИЛИЙ ВЕЛИКИЙ, ГРИГОРИЙ БОГОСЛОВ, ИОАНН ЗЛАТОУСТ, ФОМА АКВИНСКИЙ, ГРИГОРИЙ ПАЛАМА. РЕАЛИЗМ И НОМИНАЛИЗМ 45.65 KB
  которое признавало Бога в качестве высшего существующего начала а весь окружающий мир Его творения. 50 Можно выделить следующие основные черты средневековой теологической философии: теоцентризм главной причиной всего сущего высшей реальностью основным предметом философских исследований являлся Бог; изучению самого по себе космоса природы явлений окружающего мира уделялось мало внимания так как они считались творением Бога; господствовали догматы истины не нуждающиеся в доказательствах о творении всего Богом и откровении...
32877. Географическая среда в развитие общества. Географическое направление в обществоведении 41.93 KB
  Географическая среда в развитие общества. В более узком смысле – это часть материи за исключением общества т. Географическая среда имеет историческую обусловленность и расширяется с развитием общества. Географическая среда является важным фактором развития общества: Она влияет на разделение труда на размещение отраслей производства что изучается экономической географией; Она влияет на темпы развития общества может быть боле или менее благоприятной; Она опосредованно влияет на характер политических систем что изучается...
32878. Материальное производство. Труд. Производительные силы. Технологический способ труда. Технократическое направление в обществоведении 41.14 KB
  Технологический способ труда. Лишь для земледельческого труда делалось некоторое исключение. Средневековую цивилизацию многие историки культуры трактуют как цивилизацию труда. ознаменовался новым сдвигом в социальном статусе и образе труда.
32879. Производственные отношения. Социальная закономерность. Фатализм и волюнтаризм 30.27 KB
  Производственные отношения. Производственные отношения совокупность материальных и экономических отношений между людьми в процессе общественного производства и движения общественного продукта от производства до потребления. Производственные отношения являются необходимой стороной общественного производства. В процессе труда складываются отношения обусловленные потребностями технологии и организации производства например отношения между рабочими различных специальностей между организаторами и исполнителями связанные с...
32880. Общественно-экономическая формация и цивилизация 38.11 KB
  Общественноэкономическая формация в марксизме стадия общественной эволюции характеризующаяся определённой ступенью развития производительных сил общества и соответствующим этой ступени историческим типом экономических производственных отношений которые зависят от неё и определяются ею. Не существует формационных ступеней развития производительных сил которым не соответствовали бы обусловленные ими типы производственных отношений. Как ступени исторического развития человечества следующие за варварством. Как уровень ступень...
32881. Социальная структура. Классы и классовая борьба. Интеллигенция 42.81 KB
  Устойчивые связи между любыми элементами общества устойчивые взаимозависимости корреляции; 2. В широком смысле под социальной структурой понимается строение общества в целом совокупность устойчивых связей между его основными функциональными сферами экономикой политикой культурой и др. В узком смысле под социальной структурой общества понимается деление общества на различные социальные группы системы устойчивых связей между ними а также внутренняя структура различных социальных общностей. В европейской философии XVIIXVIII...
32882. Государство. Монархия и республика. Право 45.02 KB
  Государство. Государство ядром полит. Государство – это соц. Теории происхождения государство: 1.
32883. Общественное сознание и его структура. Ценности, нравственность, искусство 39.97 KB
  Общественное сознание и его структура. Общественное сознание воззрения людей в их совокупности на явления природы и социальную реальность выраженные в созданных обществом естественном или искусственном языке творениях духовной культуры социальных нормах и взглядах соц. Сознание людей может отставать от общественного бытия и не соответствовать ему. 3 Общественное сознание активно воздействует на всю жизнь общества.