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

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

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

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


 

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

42214. ИССЛЕДОВАНИЕ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ПЬЕЗОЭЛЕКТРИЧЕСКОГО ИСПОЛНИТЕЛЬНОГО УСТРОЙСТВА 1.9 MB
  Целью работы является изучение математических моделей и исследование характеристик исполнительного устройства построенного на основе пьезоэлектрического двигателя микроперемещений. Исполнительные устройства на основе пьезоэлектрических двигателей ПД позволяют получить субмикронную 107108м точность позиционирования в диапазоне перемещения до 103м и при этом обеспечить полосу пропускания свыше 1кГц. На основании приведенных выше уравнений может быть составлена структурная схема пьезоэлектрического исполнительного устройства см.
42215. ИЗУЧЕНИЕ КОНСТРУКЦИИ И ПРИНЦИПА ДЕЙСТВИЯ УГОЛЬНОГО МИКРОФОНА И ЭЛЕКТРОМАГНИТНОГО ТЕЛЕФОНА 106 KB
  Действие угольных микрофонов основано на изменении угольного порошка под влиянием звуковых колебаний воздействующих на мембрану микрофона. Устройство микрофона в упрощенном виде и способ его включения в электрическую цепь показаны на рис. Постоянная составляющая этого тока i0 является током питания микрофона; переменная составляющая разговорным током i .
42216. Огляд систем керування базами даних (СКБД) 80 KB
  Завдання Відповідно до варіанту з таблиці 1 знайти в періодичній літературі та мережі Інтернет інформацію про СКБД. У додатках наводяться формули таблиці схеми якщо вони суттєво полегшують розуміння роботи. Якщо в роботі є рисунки і таблиці які розташовані на окремих аркушах їх слід включати до загальної нумерації. Таблиці Цифровий матеріал доцільно подавати у вигляді таблиць.
42217. Нейросетевое прогнозирование. Методические указания 204 KB
  В наиболее распространенном случае ИНС обучается прогнозу на 1 отсчет времени вперед используя предыдущих значений. Другими словами на вход ИНС предъявляется вектор и требуется чтобы на выходе ИНС появилось значение: . Обучение ИНС производится по известному временному ряду .
42218. Моделирование источника заявок в системе массового обслуживания в среде Simulink 23.5 KB
  Источник генерирует последовательность однородных заявок отличающихся моментами времени появления. Интервалы времени между моментами появления заявок являются случайными величинами с известным законом распределения параметры которого остаются постоянными в течение моделируемого интервала времени . Результатом работы источника заявок является последовательность значений в пределах от нуля до .
42219. Реализация БД визуальными средствами СУБД Access 2003 358.5 KB
  В Access 2003 имеется возможность открывать таблицы, запросы, представления, сохраненные процедуры, функции и формы в режимах сводной таблицы и сводной диаграммы. Теперь анализировать данные и создавать сложные сводные таблицы и сводные диаграммы можно гораздо проще. Существует возможность сохранять представления в режимах сводной таблицы и сводной диаграммы в качестве страниц доступа к данным, которые затем может просмотреть любой
42220. Комитетные методы обучения нейронных сетей 109.5 KB
  Применение комитетных методов теоретически не хуже применения одного классификатора. Это правило часто наблюдается и на практике однако бывают случаи когда комитетная классификация работает несколько хуже одного классификатора. обучение mго классификатора зависит от результата обучения предыдущих m1 классификаторов. При этом во время обучения mго классификатора больше внимания уделяется примерам на которых чаще ошибаются предыдущие классификаторы.
42221. Решение задачи линейного программирования 146 KB
  Кабель первого типа содержит 1 телефонных b1 телеграфных и c1 фототелеграфных каналов а кабель второго типа 1 телефонных b2 телеграфных и c1 фототелеграфных каналов. Стоимость 1 км кабеля первого типа равна p1 тыс. второго типа p2 тыс. Кабель первого типа содержит 41 телефонных 3b1 телеграфных и 2c1 фототелеграфных каналов а кабель второго типа 11 телефонных 2b2 телеграфных и 5c1 фототелеграфных каналов.
42222. Аппараты защиты. Тепловое реле тока типа ТРТ-111 и автоматический выключатель типа АК-50К 53.5 KB
  Предмет исследования Объектом исследования является тепловое реле тока типа ТРТ111 и автоматический выключатель типа АК50К [1 2]. Тепловое реле тока предназначено для автоматического отключения защищаемого электротехнического объекта при наличии в цепи длительно действующих токов перегрузок. Тепловое реле тока входит составной частью в конструкцию магнитного пускателя функционально предназначенного для управления асинхронными электродвигателями. Структурная схема исполнения типовой конструкции теплового реле тока представлена на рис.