Компенсация помех, Когда остановиться? Опции

[Михайлo]

Термины:
1. Под компенсацией определённого вида помех (в идеальном случае) будем понимать вычитание из суммы полезного сигнала и помех только составляющих искомого вида помех.
2. Элемент разрешения по дальности одной частотной пачки когерентно-импульсной радиолокационной станции - это одноимённые отсчёты, полученные на выходе радиолокационного приёмника после излучения каждого из N зондирующих импульсов.

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

модуль....фаза, градусы
А(1)=10, fi(1)=5
А(2)=10, fi(2)=5
А(3)=99, fi(3)=5
А(4)=10, fi(4)=5
А(5)=10, fi(5)=5

В этом элементе разрешения находится пассивная помеха амплитудой в 10 единиц и нулевым сдвигом допплеровской фазы. Среди откликов от третьего зондирующего импульса присутствует импульсная помеха с амплитудой 89 ед. и начальной фазой в 5 градусов.

Здесь - 2 импульсные помехи (среди откликов от 3-го и 7-го зондирующих импульсов):

модуль....фаза, градусы
А(1)=10, fi(1)=150
А(2)=15, fi(2)=155
А(3)=99, fi(3)=150 ...... должно быть: А(3)=20
А(4)=25, fi(4)=165
А(5)=30, fi(5)=170
А(6)=90, fi(6)=175 ...... должно быть: А(6)=35
А(7)=40, fi(7)=180

Первый раз мне нужно вычесть из A(3) некоторым образом вычисленное значение помехи в 70 ед. и начальной фазой в 150 градусов. Второй раз - из А(6) помеху в 55 ед. и начальной фазой в 175 градусов.

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

[Михайлo]

Спасибо всем заинтересовавшимся!
Что такое медианный фильтр - доподлинно не знаю - пойду посмотрю.
Количество отсчётов функции - 8...17. Никакой доп. информации кроме входный отсчётов нет.
С аппроксимацией огибающей (по различным алгоритмам: сглаживание, линейная регрессия, сравнение очередного значения с его экстраполированным по предыдущим значениям [а вот на счёт сплайнов - уже не помню]) и последующим нахождением выпавших точек я пытался работать, но в целом удовлетворяющего меня результата не достиг.
Дело тут, наверное, вот в чём: отсчёты у меня образуют двумерную функцию. Поэтому расмотрение каждой плоскости отдельно - не лучший вариант. Сейчас я пытаюсь оценить необходимость в проведенной компенсации очередного отсчёта на выходе преобразования Фурье: пассивная помеха когерентно накопится в нулевом фильтре, а в остальных фильтрах будет присутствовать импульсная помеха (преимущественно). Но вот что-то не получается разработать универсальный критерий...
Но, может быть, кто-то предложит ещё какой-нибудь вариант?

Вот 3 примера данных. Восклицательными знаками (справа) отмечены отсчёты, в которых присутствуют импульсные помехи. Точки использованы вместо пробелов, чтобы сохранить столбцы...

A[ 0]= 35 fi= 46° !... A[ 0]= 54 fi=-121° .... A[ 0]= 47 fi= 168° !...
A[ 1]= 39 fi= 21° .... A[ 1]= 44 fi=-120° .... A[ 1]= 47 fi= 115°
A[ 2]= 56 fi= 18° !... A[ 2]= 38 fi=-117° .... A[ 2]= 89 fi= 117° !...
A[ 3]= 41 fi= 22° .... A[ 3]= 46 fi=-100° .... A[ 3]= 52 fi= 116°
A[ 4]= 40 fi= 20° .... A[ 4]= 54 fi= -97° ..... A[ 4]= 53 fi= 75° .!...
A[ 5]= 43 fi= 21° .... A[ 5]= 59 fi=-104° .... A[ 5]= 59 fi= 118°
A[ 6]= 39 fi= 24° .... A[ 6]= 59 fi=-110° .... A[ 6]= 62 fi= 118°
A[ 7]= 43 fi= 23° .... A[ 7]= 56 fi=-111° .... A[ 7]= 64 fi= 118°
A[ 8]= 43 fi= 22° .... A[ 8]= 51 fi=-107° .... A[ 8]= 68 fi= 117°
A[ 9]= 60 fi= 19° !... A[ 9]= 16 fi=-130° !... A[ 9]= 70 fi= 116°
A[10]= 45 fi= 24°.... A[10]= 49 fi=-105° .... A[10]= 73 fi= 114°
A[11]= 45 fi= 22°.... A[11]= 46 fi=-108° .... A[11]= 75 fi= 116°
A[12]= 42 fi= 21°.... A[12]= 49 fi=-104° .... A[12]= 82 fi= 117°
A[13]= 45 fi= 22°.... A[13]= 47 fi= -92° ..... A[13]= 83 fi= 114°
A[14]= 45 fi= 22°.... A[14]= 52 fi= -87° ..... A[14]= 87 fi= 115°
A[15]= 46 fi= 27°.... A[15]= 58 fi= -91° ..... A[15]= 92 fi= 116°

Сообщение отредактировал Михайлo - May 25 2007, 05:31

[-=ВН=-]

Спасибо всем заинтересовавшимся!
Что такое медианный фильтр - доподлинно не знаю - пойду посмотрю.
Количество отсчётов функции - 8...17. Никакой доп. информации кроме входный отсчётов нет.
С аппроксимацией огибающей (по различным алгоритмам: сглаживание, линейная регрессия, сравнение очередного значения с его экстраполированным по предыдущим значениям [а вот на счёт сплайнов - уже не помню]) и последующим нахождением выпавших точек я пытался работать, но в целом удовлетворяющего меня результата не достиг.
Дело тут, наверное, вот в чём: отсчёты у меня образуют двумерную функцию. Поэтому расмотрение каждой плоскости отдельно - не лучший вариант. Сейчас я пытаюсь оценить необходимость в проведенной компенсации очередного отсчёта на выходе преобразования Фурье: пассивная помеха когерентно накопится в нулевом фильтре, а в остальных фильтрах будет присутствовать импульсная помеха (преимущественно). Но вот что-то не получается разработать универсальный критерий...
Но, может быть, кто-то предложит ещё какой-нибудь вариант?

Вот 3 примера данных. Восклицательными знаками (справа) отмечены отсчёты, в которых присутствуют импульсные помехи. Точки использованы вместо пробелов, чтобы сохранить столбцы...

A[ 0]= 35 fi= 46° !... A[ 0]= 54 fi=-121° .... A[ 0]= 47 fi= 168° !...
A[ 1]= 39 fi= 21° .... A[ 1]= 44 fi=-120° .... A[ 1]= 47 fi= 115°
A[ 2]= 56 fi= 18° !... A[ 2]= 38 fi=-117° .... A[ 2]= 89 fi= 117° !...
A[ 3]= 41 fi= 22° .... A[ 3]= 46 fi=-100° .... A[ 3]= 52 fi= 116°
A[ 4]= 40 fi= 20° .... A[ 4]= 54 fi= -97° ..... A[ 4]= 53 fi= 75° .!...
A[ 5]= 43 fi= 21° .... A[ 5]= 59 fi=-104° .... A[ 5]= 59 fi= 118°
A[ 6]= 39 fi= 24° .... A[ 6]= 59 fi=-110° .... A[ 6]= 62 fi= 118°
A[ 7]= 43 fi= 23° .... A[ 7]= 56 fi=-111° .... A[ 7]= 64 fi= 118°
A[ 8]= 43 fi= 22° .... A[ 8]= 51 fi=-107° .... A[ 8]= 68 fi= 117°
A[ 9]= 60 fi= 19° !... A[ 9]= 16 fi=-130° !... A[ 9]= 70 fi= 116°
A[10]= 45 fi= 24°.... A[10]= 49 fi=-105° .... A[10]= 73 fi= 114°
A[11]= 45 fi= 22°.... A[11]= 46 fi=-108° .... A[11]= 75 fi= 116°
A[12]= 42 fi= 21°.... A[12]= 49 fi=-104° .... A[12]= 82 fi= 117°
A[13]= 45 fi= 22°.... A[13]= 47 fi= -92° ..... A[13]= 83 fi= 114°
A[14]= 45 fi= 22°.... A[14]= 52 fi= -87° ..... A[14]= 87 fi= 115°
A[15]= 46 fi= 27°.... A[15]= 58 fi= -91° ..... A[15]= 92 fi= 116°

Как вариант - привлечь кепстры Спектр мощности от логарифма спектра мощности. В кепстре все импульсные помехи, вместе с белым шумом, сосредоточатся в области нулевых сачтот. Можно там их и похерить лифтром соответствующим.
Еще вариант, наверное уже упоминавшийся, аппроксимация данных, или возможно их куска, какой-то простой функцией, типа полинома невысокой степени. Возможно раздельно для амлитуды и фазы. Вычитание из данных найденной аппроксимирующей функции. В полученной разности (у нее видимо будет в основном помехо-шумовой характер) - определение положения импульсных помех и их последующее убиение в исходных данных.
Определения положения можно так попробовать.
Оценить дисперсию шума. В самом простом случае - прямо по полученной разности.
А именно - найти медиану квадрата разности, медиана даст в данном случае нечто близкое ко 2 моменту шума. А шум скорее всего будет с близким к 0 матожиданием, поэтому 2 момент - фактически дисперсия. По дисперсии поставить порог. Для чего умножить ее на коэффициент. Коэффициент - из вероятноти тревоги. Оч. часто к-т ~=20-25. Результат в качестве порога. Ищете максимум в квадрате разности. Сравниваете с порогом . Превышает - помеха. Запоминаете ее положение и величину и убиваете в квадрате разности. Опять ищете максмум и делаете то же самое. Так до тех пор, пока превышения порога не прекратятся. Затем корректируете исх. данные по найденным положениям и величинам.
P.S. суммарная длительность (или в Вашем случае скорее число) помех должна быть меньше не менее чем в 2 раза длительности (размера) окна анализа.

Сообщение отредактировал -=ВН=- - May 25 2007, 15:01

[Stanislav]

Как вариант - привлечь кепстры Спектр мощности от логарифма спектра мощности. В кепстре все импульсные помехи, вместе с белым шумом, сосредоточатся в области нулевых сачтот. Можно там их и похерить лифтром соответствующим...

Простите, не понял. Дальше-то что с ним делать? Если не ошибаюсь, автору нужно восстановить именно форму зашумлённых функций...

Вот, программку набросал на матлабе. Можно "не отходя от кассы" проверить её действие.

Код

function puls_remove()

% Задаём сетки
f_len=16;
sampl = [1:f_len]';
sampl1=sampl;

% Задаём порог
Thresh=20;
Thresh1=2.5;

% Задаём значения функций
A = [35 39 56 41 40 43 39 43 43 60 45 45 42 45 45 46]';
fi= [46 21 18 22 20 21 24 23 22 19 24 22 21 22 22 27]';

% A = [54 44 38 46 54 59 59 56 51 16 49 46 49 47 52 58]';
% fi= [-121 -120 -117 -100 -97 -104 -110 -111 -107 -130 -105 -108 -104 -92 -87 -91]';

% A = [47 47 89 52 53 59 62 64 68 70 73 75 82 83 87 92]';
% fi= [168 115 117 116 75 118 118 118 117 116 114 116 117 114 115 116]';

% Аппроксимируем функции полиномами первого порядка (отрезками прямых)
fit_A = fit(sampl, A, 'poly1');
A_fit = feval(fit_A, sampl);
fit_fi= fit(sampl, fi,'poly1');
fi_fit = feval(fit_fi, sampl);

% Рисуем функции и аппроксимации
figure(1)
plot (A, '.-', 'Color', [0 0.6 0], 'MarkerEdgeColor','k')
hold on
plot (A_fit, '.-', 'Color', [0.9 0 0.9], 'MarkerEdgeColor','k')
grid on
set(gca,'XTick',1:16)

figure(2)
plot (fi, '.-', 'Color', [0 0.6 0], 'MarkerEdgeColor','k')
hold on
plot (fi_fit, '.-', 'Color', [0.9 0 0.9], 'MarkerEdgeColor','k')
% hold off
grid on
set(gca,'XTick',1:16)


% Вычисляем модули расстояний
dist_A =abs(A-A_fit);
dist_fi=abs(fi-fi_fit);

%********************************************************************
% Вычисляем их медианы
med_A = median(dist_A);
med_fi= median(dist_fi);

% Нормализуем (приводим медианы к одному и тому же значению)
dist_fi = dist_fi * med_A/med_fi;

% Вычисляем квадрат эвклидова расстояния в каждой точке выборки
dist2 = dist_A.*dist_A + dist_fi.*dist_fi;

% Отбрасываем ненужные точки
for (i=16:-1:1)
    if (dist2(i) > Thresh*med_A*med_A)
        A(i)=[]; fi(i)=[]; sampl1(i)=[];
    end
end

%*********************************************************************
% % Или вычисляем среднюю мощность (требует другого порога!)
% msq_A = (dist_A'*dist_A)/f_len;
% msq_fi = (dist_fi'*dist_fi)/f_len;
%
% % Нормализуем (приводим СКЗ к одной и той же величине)
% dist_fi = dist_fi * sqrt(msq_A/msq_fi);
%
% % Вычисляем квадрат эвклидова расстояния в каждой точке выборки
% dist2 = dist_A.*dist_A + dist_fi.*dist_fi;
%
% % Отбрасываем ненужные точки
% for (i=16:-1:1)
%     if (dist2(i) > Thresh1*msq_A)
%         A(i)=[]; fi(i)=[]; sampl1(i)=[];
%     end
% end
%*********************************************************************

% Делаем аппроксимацию функций с выброшенными точками полиномами
% 2-го (или иного) порядка.
fit_A = fit(sampl1, A, 'poly2');
A_fit = feval(fit_A, sampl);
fit_fi= fit(sampl1, fi,'poly2');
fi_fit = feval(fit_fi, sampl);

% Рисуем окончательный результат (красным-функции с выброшенными точками)
figure(1)
plot (sampl1, A, '.-', 'Color', [0.8 0 0], 'MarkerEdgeColor','k')
plot (sampl, A_fit, '.-', 'Color', [0 0 1], 'MarkerEdgeColor','k')
hold off

figure(2)
plot (sampl1, fi, '.-', 'Color', [0.8 0 0], 'MarkerEdgeColor','k')
plot (sampl, fi_fit, '.-', 'Color', [0 0 1], 'MarkerEdgeColor','k')
hold off

%*********************************************************************
% Окончательные значения полиномов A_fit и fi_fit можно взять в качестве
% выходных значений функций (синие линии).
%*********************************************************************
return;

Прогга - именно набросок, а не истина в последней инстанции, и является лишь демонстрацией принципа.

PS. Вместо средней мощности лучше взять её медианное значение, как и предложил уважаемый BH.

[-=ВН=-]

Простите, не понял. Дальше-то что с ним делать? Если не ошибаюсь, автору нужно восстановить именно форму зашумлённых функций...

А дальше обратно все проделать (после вырезания околонулевой области). Фаза предварительно сохранена.
Вариант дан больше шутки ради, он трудоемок вычислительно. Но вообще он рабочий, хотя естественно полной очистки не дает, да и тонкостей там хватает. Когда-то давно я его пробовал.
Как раз для борьбы с широкополосными помехами.