Главная страница | Контакты | Случайная страница Автомобили Астрономия Биология География Дом и сад Другие языки Другое Информатика История Культура Литература Логика Математика Медицина Металлургия Механика Образование Охрана труда Педагогика Политика Право Психология Религия Риторика Социология Спорт Строительство Технология Туризм Физика Философия Финансы Химия Черчение Экология Экономика Электроника

Расчёт отношения сигнал/шум

(где t – время наличия полезного сигнала)

дБ

ОСШ в моём случае 17 дБ (и колеблется около этого значения при многократных опытах)

Код программы (с достаточно подробными комментариями)

clear all;

close all;

f_rel = 0.005; % 200 отсчётов на период синуса

T = 1000; % Время моделирования

S_1(1:T) = 10 * cos (f_rel * 2 * pi * (0:T-1)); % сигнал из 1го канала

S_1_hilb(1:T) = 10 * sin (f_rel * 2 * pi * (0:T-1)); %сигнал из 1го канал преобр. по Гильберту

S_2(1:T) = 10 * cos (f_rel * 2 * pi * (100:T+99)); % сигнал из 2го канала

% Предположим, что сигнал есть в периоды от T = 100..250 и 650..800

ff1_1 = 1; ff2_1 = 251; ff3_1 = 801;

ff1_2 = 99; ff2_2 = 549; ff3_2 = 1000;

true_impulse(1:T) = 1;

S_1(ff1_1:ff1_2) = 0; S_2(ff1_1:ff1_2) = 0; true_impulse(ff1_1:ff1_2) = 0;

S_1(ff2_1:ff2_2) = 0; S_2(ff2_1:ff2_2) = 0; true_impulse(ff2_1:ff2_2) = 0;

S_1(ff3_1:ff3_2) = 0; S_2(ff3_1:ff3_2) = 0; true_impulse(ff3_1:ff3_2) = 0;

% создаём массив белого Гауссовского шума мощностью 0дБ:

noise_1 = wgn(T,1,0);

noise_2 = wgn(T,1,0);

% чтобы сделать зашумлённый сигнал нужно сложить массивы одинаковых размерностей, для этого транспонируем вектор столбец

noise_1 = noise_1';

noise_2 = noise_2';

% Создаём зашумлённые сигналы:

S_noise_1 = S_1 + noise_1;

S_noise_hilb_1 = S_1_hilb + noise_1; % принципиально ли, что шум не проходит через преобразователь Гильберта (в модели конечно же)???

S_noise_2 = S_2 + noise_2;

%------------- конец формирования сигналов ------------

%------------- выделение энергии сигнала --------------

X_Energy_noise = S_noise_1.* S_noise_2;

Y_Energy_noise = S_noise_hilb_1.* S_noise_2;

% энергия сигнала в текущий момент времени (или шума, если сигнала нет):

A = sqrt(X_Energy_noise.^2 + Y_Energy_noise.^2);

%---------------объявление переменных ---------------------

energy_stack(1:11) = 0; % стек, если 9 или 10 единиц - сигнал есть

threshold = 10; % порог срабатывания на энергию сигнала

detection_impulse(1:T) = 0; %импульс, говорящий о наличии сигнала

%--------------- переменные, относящиеся к "двойному условию" ------

unstoppable_counter = 0;

% возможность фронта при определённом количестве непрерывных единиц:

probauble_front_number = 4;

label(1:T) = 0; % метка о возможном фронте

%---------------конец объявления переменных -------------------------

%-----заполнения стека по энергии отсчёта и синтез импулься обнаружения ---

for i=3:T

%занесение единицы в стек если энергия больше порога

if A(i)>threshold

energy_stack(1) = 1;

else energy_stack(1) = 0;

end

%определение наличия 9 или 10 единиц в стеке

counter = 0; % Обнуление счётчикаж

for n = 1:10

if energy_stack(n) == 1

counter = counter + 1;

end

end

% если 9 или 10 единиц - создаётся импульс обнаружения

if counter > 8

detection_impulse(i) = 1;

end

%----------------------------------------------------------------

%- возможное досрочное обнаружение сигнала при помощи двойного условия

% мощности и меньшего количества единиц.в стеке энерг. обнаружителя

if A(i) > threshold

unstoppable_counter = unstoppable_counter + 1;

else unstoppable_counter = 0;

end

if unstoppable_counter > probauble_front_number

label(i) = 1;

end

%------------------------------------------------------------------

%сдвиг стека вправо >>

for n = 1:9

energy_stack(11-n) = energy_stack(10-n);

end

% при кратковременном прерывании импульса считается, что импульс есть и

% он не прерывался

if detection_impulse(i-1) == 0

if (detection_impulse(i-2) == 1) && (detection_impulse(i) == 1);

detection_impulse(i-1) = 1;

end

end

end

%-------------- "Мощностная часть" --------------------

%-------------объявление переменных -------------------

P(1:T) = 0; % мгновенная мощность сигнала/шума

m = 4; % количество ячеек в стеке мощности

power_stack(1:m) = 0; % стек, говорящий об изменении мощности

power_threshold = 5;

power_det_imp(1:T) = 0; % импульс, говорящий о изменении мощности

% расчёт мощности

for i=2:T

P(i) = (A(i) - A(i-1));

end

test(1:T) = 0;

for i=18:T

P_round(i) =abs((P(i) + P(i-2) + P(i-4) + P(i-6) + P(i-8) + P(i-10) + P(i-12) + P(i-14) + P(i-16))/9); % поиск производной по 18 точкам

%занесение в стек 1, если изменение мощности больше 1

if P_round(i) > 0.9 * P_round(i-m)

power_stack(1) = 1;

else power_stack(1) = 0;

end

%определение кол-ва единиц в стеке:

counter = 0;

for n = 1:m

if power_stack(n) == 1

counter = counter + 1;

test(i) = counter;

end

end

%если m-1 единиц в стеке - сигнал есть

if (counter > m-2) && (label(i)==1)

power_det_imp(i) = 1;

end

%сдвиг стека вправо >>

for n = 1:(m-1)

power_stack((m+1)-n) = power_stack(m-n);

end

end

%--- объединяя условия наличия сигнала по мощности и энергии

double_det_imp(1:T) = 0;

for i=1:T

if power_det_imp(i)==1 || detection_impulse(i)==1

double_det_imp(i) = 1;

else double_det_imp(i) = 0;

end

end

%--------------------------------------------------

%--------------- расчёт отношения Сигнал - Шум.. -------------

noise_sum = 0; S_sum = 0;

noise_sq = 0; S_sq = 0;

% расчёт времени полезного сигнала (чтобы знать, на что делить при расчёте.

signal_time = T - ((ff1_2 - ff1_1 + 1) + (ff2_2 - ff2_1 + 1) + (ff3_2 - ff3_1 + 1));

for i=1:T %нахождение суммы квадратов сигнала и шума.

noise_sum = noise_1(i)^2 + noise_2(i)^2 + noise_sum;

S_sum = S_1(i)^2 + S_sum;

end

noise_sq = sqrt (noise_sum / (2 * T)); % среднеквадратическое знач. шума

S_sq = sqrt (S_sum / signal_time); % Среднеквадр. значение сигнала

OSW = 20 * log10(S_sq/noise_sq); % вычисление ОСШ

%----------------- построение графиков ---------------------------------

figure(1);

plot(1:T, S_noise_1, 1:T, 10 * true_impulse, 1:T, 8 * detection_impulse);

legend('Исходный сигнал','априорный импульс наличия сигнала', 'импульс обнаружения с энергетического обнаружителя');

title(' график 1 ');

xlabel(' отсчёты i ');

ylabel(' U, mV ');

figure(2);

plot(1:T, P_round, 1:T, 10 * power_det_imp, 1:T, 20 * true_impulse);

legend('Мощность, вычесленная по 18 отсчётам', 'импульс обнаружения с помощью мощности и энергии', 'априорный импульс наличия сигнала');

title('график 2 ');

xlabel(' отсчёты i ');

ylabel(' P, mW ');

figure(3);

plot(1:T, 3 * true_impulse, 1:T, 1 * detection_impulse, 1:T, 2 * double_det_imp);

legend(' априорный импульс наличия сигнала ','импульс обнаружения с энергетического обнаружителя', '"сдвоенный по функции или" импульс обнаружения');

title(' график 3 ');

xlabel(' отсчёты i ');

ylabel(' U, V ');