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

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

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

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


 

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

74838. Житие Сергия Радонежского. Общая характеристика агиографического творчества Пахомия Лагофета 16.33 KB
  Епифаний создал Житие Сергия Радонежского. Епифаний хорошо передает факты биографии Сергия с лирической теплотой говорит о его деятельности связанной с борьбой против ненавистной розни за укрепление централизованного Русского государства. О роли Сергия Радонежского и Стефана Пермского в политическом и нравственном возрождении Русской земли говорил В.
74839. Повесть о взятии Царьграда Нестора – Искандера. Историко-философское осмысление событий в повести 15.23 KB
  Повесть содержит описание истории Константинополя с момента его основания, но особенно подробно рассказывается об осаде византийской столицы турками и ее взятии. Хотя в “Повести” содержится немало достоверных сведений, в целом она все же чисто литературное произведение, а не документальная хроника. Некоторые эффектные сюжетные коллизии оказываются вымыслом: так, в Константинополе во время осады не было патриарха
74840. Хождение за три моря тверского купца Афанасия Никитина 17.94 KB
  Хождение за три моря Афанасия Никитина. является Хождение за три моря тверского купца Афанасия Никитина помещенное под 1475 г. Да станет Русская земля благоустроенной и да будет в ней справедливость Православная вера является для Никитина символом родины.
74841. Памятник Муромо-Рязанской литературы – «Повесть о Петре и Февронии» Поэтичность, демократизм повести. Образ крестьянки Февронии 17.58 KB
  Герои повести — исторические лица: Петр и Феврония княжили в Муроме в начале XIII века, они умерли в 1228 г. Однако в повести историчны только имена, вокруг которых был создан ряд народных легенд, составивших основу сюжета повести. Как указывает М. О. Скрипиль, в повести объединены два народнопоэтических сюжета: волшебной сказки об огненном змее и сказки о мудрой деве.
74842. Домострой – книга, утвердившая нормы семейной жизни 15.89 KB
  Следует отнести также Домострой составление которого приписывается благовещенскому попу Сильвестру входившему в Избранную раду. Домострой регламентировал поведение человека как в государственной так и в семейной жизни.
74843. Бытовая повесть. Повесть о Горе и Злочастии 17.34 KB
  Основу сюжета повести составляет трагическая история жизни Молодца отвергнувшего родительские наставления и пожелавшего жить по своей воле как ему любо. Постоянная опека родителей не научила Молодца разбираться в людях понимать жизнь и он платится за свою доверчивость за слепую веру в святость уз дружбы. Губит Молодца похвальба своим счастьем и богатством. В сознании Молодца еще живучи традиционные представления.
74844. Житие протопопа Аввакума, им самим написанное. Автобиография «Жития». Демократизм и публицистичность «Жития». Традиция и новаторство в жанре Жития 21.16 KB
  Житие протопопа Аввакума им самим написанное. Аввакум так определяет рамки своего повествования. Центральная тема жития тема личной жизни Аввакума неотделимая от борьбы за древлее благочестие против Никоновых новшеств.
74845. Публицистика петровского времени. Юности честное зерцало 14.87 KB
  Зерцало было издано в соответствии с духом петровских реформ когда основу всей книгопечатной продукции составляли разного рода руководства и наставления. Вторая часть это собственно зерцало то есть правила поведения для младых отроков и девушек дворянского сословия.
74846. Екатерина II. «Всякая всячина». Н.И. Новиков, его литературные труды и сатирические журналы. Сатира на лица и на пороки 16.52 KB
  Новиков его литературные труды и сатирические журналы. Предполагают что свои материалы в журнал посылали Фонвизин князь Щербатов и даже Новиков. Новикова. В Трутне под подписью Правдолюбов Новиков упрекает Екатерину.