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

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

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

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


 

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

25732. Спутниковые системы связи. Классификация ИСЗ по особенностям орбиты. Спутниковые службы в системах связи 15.05 KB
  Классификация ИСЗ по особенностям орбиты. Использование ИСЗ позволяет резко повысить дальность радиосвязи тк ретранслятор располагается высоко над Землей. 3 основных вида ИСЗ: ИСЗ на высокой эллиптической орбите ВЭО ИСХ на геостационарной орбите ГЭО ИСЗ на низковысотной орбите НВО ВЭО Спутники типа молния с периодом обращения 12 часов наклоном орбиты 63 градуса высотой апогея над северным полушарием 40 тыс. В области апогея скорость движения ИСЗ замедляется и обеспечивается радиовидимость 68 часов.
25733. Распространение декаметровых волн 37.72 KB
  К диапазону KB декаметровые волны относят радиоволны с длиной волны от 100 до 10м. В отличие от более коротких волн которые распространяются земной волной декаметровые волны распространяются в основном путем отражении от ионосферы. Радиус действия земной волны в диапазоне коротких волн сравнительно невелик и при обычно используемых мощностях передатчиков не превышает нескольких десятков километров. Но короткие волны могут распространяться на многие тысячи километров путем многократных последовательных отражений от ионосферы и Земли и...
25734. Взаимодействие уровней модели OSI 23.42 KB
  Каждый уровень модели OSI выполняет определенную задачу в процессе передачи данных по сети. Уровень 7 Applicayion layer A Прикладной Ур. Каждый уровень компьютераотправителя взаимодействует с таким же уровнем компьютераполучателя как будто он связан напрямую. Каждый уровень модели выполняет свою функцию.
25735. Каналы связи. Классификация каналов связи. Параметры каналов связи. Условие передачи сигнала по каналу связи 287 KB
  Канал связи — система технических средств и среда распространения сигналов для передачи сообщений (не только данных) от источника к получателю (и наоборот). Канал связи, понимаемый в узком смысле (тракт связи), представляет только физическую среду распространения сигналов, например, физическую линию связи.
25737. Рынок современной прессы 34.17 KB
  газет порядка 18 тыс. Преобладающим товаром и основным информационным носителем на рынке российской периодики как и прежде остаются газеты общий тираж которых в 2005 году составил 8 млрд 312 млн экз. из которых 29 млрд приходится на долю общероссийских газет. В отличие от газет подавляющее большинство журналов 916 наименований и 98 совокупного годового тиража приходится на Москву и СанктПетербург.
25738. Законодательство в сфере журналистики. Закон о СМИ 25.88 KB
  Право на получение информации одновременно означает обязанность государственных органов и органов местного самоуправления всех властных и общественных структур дать ответ на обращение к ним гражданина. Право производить и распространять информацию представляет собой творческий процесс по созданию информации в любой форме от книги до сигналов в космические дали и свободному ее распространению в том числе с помощью различных технических средств. Недопустимость идеологической ангажированности средств массовой информации закреплена в статье...
25739. Журналистское произведение как текст – жанровые разновидности 29.76 KB
  В прессе предлагают разделять тексты на 5 групп: оперативноновостные заметка оперативноисследовательские интервью репортажи отчеты исследовательсконовостные комментарий рецензия корреспонденция исследовательские статья письмо обозрение исследовательскообразные очерк фельетон. Информационные новость заметка интервью репортаж отчет реплика Аналитические корреспонденция статья комментарии рецензия обозрение Художественнопублицистические очерк зарисовка эссе сатирические жанры – памфлет...
25740. Явления демассификации и дигитализации СМИ 19.61 KB
  Явления демассификации и дигитализации СМИ Под воздействием новых информационных технологий происходят значительные изменения в области массовых коммуникаций. Рассмотрим основные процессы лежащие в основе модификации современной системы СМИ. Под этим подразумевается перевод содержания СМИ во всех его формах – текстовой графической звуковой – в цифровой формат понятный современным компьютерам. Дигитализация устраняя различия между отдельными СМИ и уравнивая их содержание прокладывает дорогу их конвергенции.