Решить задачу обнаружения сигнала

25.11.2020

Уважаемые представители сообщества! Не получаеется постороить согласованный фильтр, АКФ функций для временной фильтрации и обратное преобразование Фурье, график. В файле лабораторная работа необходимо проделать п2.-п5 в среде matlab, а не в маткаде, не могу понять в чем проблема... Также прикладываю скрины работы в маткаде... Я буду очень признателен за помощь!

Код на матлабе:

% % узкополосная фильтрация

x0 = 1.5;
w0 = 2*pi*60;
t = [0:0.001:0.8]';
N = 4096;
 
n = length(t);
Y = zeros(1,n);
for i = 1:n
    if t(i) >= 0 && t(i) <= 0.4
        Y(i) = 0;
    elseif t(i) > 0.4  && t(i) <= 0.8 
        Y(i) = x0*cos(w0*t(i));
    end
end
plot (t,Y)
 
% проводим дискретизацию по времени
n1 = [0:N-1];
tt = n1.*(t(end)/N);
for i = 1:length(n1)
    if tt(i) >= 0 && tt(i) <= 0.4
        Y1(i) = 0;
    elseif tt(i) > 0.4  && tt(i) <= 0.8 
        Y1(i) = x0*cos(w0*tt(i));
    end
end
figure,plot(tt,Y1)
 
% накладываем гауссовский белый шум
sigma = 2;
Z = Y1+normrnd(0,sigma,size(n1));
figure,plot(tt,Z,'k-')
 
% подвергаем преобразованию Фурье
%L = [1:N-1]';
%nfft = length(L);
ns = (length(tt));
ks=1/tt(2);
YY = (abs(fft(Y1))/(ns));
YY=YY(1:ns/2+1);
YY(2:end-1) = 2*YY(2:end-1);
ZZ = (abs(fft(Z))/(ns));
ZZ=ZZ(1:ns/2+1);
ZZ(2:end-1) = 2*ZZ(2:end-1);
k=ks*(0:ns/2)/ns;
figure,
grid on;
plot(k,YY/2);
hold on
stem(k,ZZ/2);
hold off;
xlim([0 120])

% обнуляем те элменты к спектральной функции...
%  к, где спектр полезного сигнала  не  сильно отличается от нуля, и
%  совершаем обратное преобразование Фурье
YY1 = zeros(1,length(ns));
for i = 1:length(ns)
    if ns(i) >= 40 && ns(i) <= 80
        YY1(i) = YY;
    elseif ns(i) < 40  && ns(i) > 80 
        YY1(i) = 0;
    end
end
% обатное преобразование
y = ifft(YY,ns);
figure
plot(tt,Z,'k-',tt,y+6,'r--')  %,tt,y+6,'r--')


% строим две АКФ
ZZ1 = Z-Y1;

RR = zeros(1,(ns-1));
for i = 1:ns
    RR(i) = Z(i)*Z(i+1);
end
RR = (1./(N-i))*RR;


RR1 = zeros(1,(ns-1));
for i = 1:ns-1
    RR1(i) = ZZ1(i)*ZZ1(i+1);
end
RR1 = (1./(N-i))*RR1;

figure
plot (i,RR,i,RR1);


clear
clc
TMAX=1;
N=4096;
dt=TMAX/N;
tn=0:dt:TMAX;
A=1;
% Сигнал
% sig=A*(tn>= 0.5)-(A*(tn >= 1));  %+(-A*(tn >= 1));
sig=A*(tn>=0.5).*(exp(-(5*(tn-0.5).*(tn>=0.5))));
% Сигнал + помеха
x=sig+1*randn(N+1,1)'; % сигнал с шумом
x1 = x-sig; % шум без полезного сигнала
L=length(sig); %длина паттерна
N=length(x); %длина входного сигнала
y=zeros(1,N);
y1=zeros(1,N);
for i=1:N
    for p=1:L
        if (i-p)>0
        %y(i)=y(i)+x(i-p)*sig(N-p)*dt;
%         y(i)=y(i-p)+x(p)*dt;
%         y1(i)=y(i-p)+x1(p)*dt;
        y(p)=y(i-p)+x(p)*dt;
        y1(p)=y(i-p)+x1(p)*dt;
        end
    end
end
% вывод результата
subplot (4,1,1); plot(tn,sig);grid;title('signal') % шаблон сигнала
subplot (4,1,2); plot(tn,x);grid;title('signal+noise') % реальный сигнал
subplot (4,1,3); plot(tn,y);grid;title('filter output') % выход фильтр
subplot (4,1,4); plot(tn,y1);grid;title('filter output without noise') % выход фильтра без учета шума

Лучший ответ

aBoomest+942.89

8.12.2020 12:00

Ответы

aBoomest+942.89

26.11.2020 05:40

1. А картинки в маткаде априори правильные?
2. До БПФ поправил. Дальше смотрите, правильнные картинки или как? Если нет, то желательно конкретные вопросы.
- Н/Д0.00
  
  26.11.2020 05:51
  
  То что в маткаде, то да, это правильные картинки, ошибки у меня с корреляционной функцией, обратным БПФ и согласованным фильтром....
- Н/Д0.00
  
  26.11.2020 05:55
  
  Просто когда я пытаюсь совершить ОБПФ над сигналом, у которого обнулены уже спектральные составляющие, слабо отличающиеся от нуля, то построить его он не может...
aBoomest+942.89

26.11.2020 09:07

Не, все должно работать. Скольки процентные (отн. макс.) значения вы обнуляете?
- Н/Д0.00
  
  26.11.2020 09:37
  
  ну если судить по графику прямого БПФ, то обнуляю с 43 по 76 элемент k
Н/Д0.00

26.11.2020 09:36
В методичке там есть пункты, распишу подробнее, чтобы было понятно что и зачем:

пункт 2.1 - 2.2. Задаем сигнал, проводим его дискретизацию, для выполнения дальнейших преобразований (здесь проблем нет)
```
clc; 
close all; 
clear all;
% узкополосная фильтрация
x0 = 1.5;
w0 = 2*pi*60;
t = transpose(0:0.001:0.8);
N = 4096;
 
n = length(t);
Y = zeros(n,1);
for ii = 1:n
    if t(ii) >= 0 && t(ii) <= 0.4
        Y(ii) = 0;
    elseif t(ii) > 0.4  && t(ii) <= 0.8 
        Y(ii) = x0*cos(w0*t(ii));
    end
end
plot (t,Y)

% проводим дискретизацию по времени
n1 = transpose(0:N-1);
tt = n1.*(t(end)/N);
Y1 = zeros(n,1);
for ii = 1:length(n1)
    if tt(ii) >= 0 && tt(ii) <= 0.4
        Y1(ii) = 0;
    elseif tt(ii) > 0.4  && tt(ii) <= 0.8 
        Y1(ii) = x0*cos(w0*tt(ii));
    end
end
figure,plot(tt,Y1)
```
пункты 2.3 - 2.4 добавляем к сигналу гауссовский шум, и потом производим преобразование Фурье, здесь тоже все более-менее правильно...
```
% накладываем гауссовский белый шум
sigma = 2;
sz = size(n1);
Noize = normrnd(0,sigma,sz(1),sz(2));
Z = Y1+Noize;
figure,plot(tt,Z,'k-')
 
 
% подвергаем преобразованию Фурье
%L = [1:N-1]';
%nfft = length(L);
ns = (length(tt));
ks=1/tt(2);
YY = (abs(fft(Y1))/(ns));
YY=YY(1:ns/2+1);
YY(2:end-1) = 2*YY(2:end-1);
ZZ = (abs(fft(Z))/(ns));
ZZ=ZZ(1:ns/2+1);
ZZ(2:end-1) = 2*ZZ(2:end-1);
k=ks*(0:ns/2)/ns;
figure,
grid on;
plot(k,YY/2);
hold on
stem(k,ZZ/2);
hold off;
xlim([0 120])
```
Затем, необходимо изучить, на каких частотах или при каких значениях индекса амплитудный спектр полезного сигнала существенно отличен от нуля, и обнулить спектральную функцию смеси сигнала с шумом на всех остальных частотах; выполнить обратное быстрое преобразование Фурье к полученной
таким образом спектральной функции. По итогу, я должен получить картинку, как в 3.png... Здесь ошибка в том, что график получается не совсем тот, что надо:

Т.е не получается сделать так, как в маткаде.

Затем идет пункт 3, посвященный врменной фильтрации в блоке АКФ не получается построить функции, так как на 4.png... Вот код этого фрагмента:
```
% строим две АКФ
ZZ1 = Z-Y1; % чистый шум

kk = 1:ns/2;
RR = zeros(1,(ns/2));
for i = 1:ns/2
    RR(i) = Z(i)*Z(i+1);
end
RR1 = (1./(N-i))*RR;


RR1 = zeros(1,(ns/2));
for i = 1:ns/2
    RR1(i) = ZZ1(i)*ZZ1(i+1);
end
RR11 = (1./(N-i))*RR1;

figure
plot (kk,RR1,kk,RR11);
xlim([0 500]);
```
В согласованной фильтрации, для импульса, который начинается на интервале 0,5с (середина интервала) выход фильтра должен иметь такую форму (10.png), или же я что-то не понимаю...

Для эксопоненциального импулсьа это должно выглядить так (11.png)

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

Задайте смесь полезного сигнала и шума аналогично вышеизложенному.

Импульсную характеристику согласованного фильтра постройте
по образцу полезного сигнала: w(t)=x(T-t). Импульсную характеристику также подвергните дискретизации.

Постройте свёртку смеси сигнала с шумом и только с шумом. Сравните получаемые значения. Можно ли по свёртке увидеть, в какой момент начинается полезный сигнал?

Вот код этого фрагмента
```
clear
clc
TMAX=1;
N=4096;
dt=TMAX/N;
tn=0:dt:TMAX;
A=1;
% Сигнал
% sig=A*(tn>= 0.5)-(A*(tn >= 1));  %+(-A*(tn >= 1));
sig=A*(tn>=0.5).*(exp(-(5*(tn-0.5).*(tn>=0.5))));
% Сигнал + помеха
x=sig+1*randn(N+1,1)'; % сигнал с шумом
x1 = x-sig; % шум без полезного сигнала
L=length(sig); %длина паттерна
N=length(x); %длина входного сигнала
y=zeros(1,N);
y1=zeros(1,N);
for i=1:N
    for p=1:L
        if (i-p)>0
        %y(i)=y(i)+x(i-p)*sig(N-p)*dt;
%         y(i)=y(i-p)+x(p)*dt;
%         y1(i)=y(i-p)+x1(p)*dt;
        y(p)=y(i-p)+x(p)*dt;
        y1(p)=y(i-p)+x1(p)*dt;
        end
    end
end
% вывод результата
subplot (4,1,1); plot(tn,sig);grid;title('signal') % шаблон сигнала
subplot (4,1,2); plot(tn,x);grid;title('signal+noise') % реальный сигнал
subplot (4,1,3); plot(tn,y);grid;title('filter output') % выход фильтр
subplot (4,1,4); plot(tn,y1);grid;title('filter output without noise') % выход фильтра без учета шума
```
Н/Д0.00

26.11.2020 09:36
aBoomest+942.89

29.11.2020 19:55
```
YY1 = zeros(1,length(ns));
for i = 1:length(ns)
    if ns(i) >= 40 && ns(i) <= 80
        YY1(i) = YY;
    elseif ns(i) < 40  && ns(i) > 80 
        YY1(i) = 0;
    end
end
```
Вы вот тут делаете обрезание не значимых частот, если я правильно понял.
```
% обатное преобразование
y = ifft(YY,ns);
figure
plot(tt,Z,'k-',tt,y+6,'r--')  %,tt,y+6,'r--')
```
Далее ОБПФ от спектра YY.
```
YY = (abs(fft(Y1))/(ns));
```
Однако выше YY - это у вас почему то модуль спектра, а не сам спектр. Вы модуль берите потом когда рисуете график, либо промежуточную переменную дополнительную надо. А сам спектр - это совсем не модуль.
- Н/Д0.00
  
  29.11.2020 20:24
  Надо построить именно модуль преобразования, я тогда модуль возьму, когда график буду строить, а что насчёт построение двух АКФ? Вот не получается переписать маткад языком матлаба, и попытался применить функцию autocorr, графики получились похожи, но по амплитуде, немного не соответствуют тому, что получилось в маткад... Как решить эту проблему?
  
  По-поводу преобразования, я Вас понял! Я приложу полный код в виде m.файла, там и через autocorr, и через прямую запись, через сумму... а насчёт преобразования, я тогда исправлю...
```
Corr_auto1 = autocorr(Z);
% Corr_auto2 = autocorr(Z1);
% figure (6)
% r1 = autocorr(Z,'NumLags',300,'NumSTD',2);
% r2 = autocorr(Z1,'NumLags',300,'NumSTD',2);
% plot (1:length(r1),r1,1:length(r2),r2)
% ylim([-1 1]);

% plot (length(Corr_auto1),Corr_auto1,length(Corr_auto2),Corr_auto2);
```
- Н/Д0.00
  
  29.11.2020 20:31
  
  Ой, но так тогда у меня не получается переписать, чтобы сначала был спектр, затем нанести его модуль на график, а потом уже вычленить незначащие частоты и построить результат ОБПФ вместе с исходным зашумленным сигналом, т.е. как на скринах в маткаде
  - aBoomest+942.89
    
    30.11.2020 10:35
    
    Немного не так.
    
    1. У вас д.б. спектр - как есть. Спектр - есть спектр.
    
    2. По амплитуде спектра вы должны определить элементы массива подлежащие "занулению".
    
    3. И занулить их надо именно в спектре, а не в модуле, иначе от чего делать обратную процедуру.
    А отрисовывать конечно надо либо амплитуду и фазу, либо re и im, это уж как удобно.
  - Н/Д0.00
    
    30.11.2020 10:52
    
    Попытаюсь конечно переписать, хотя 3 недели уже никак не могу сделать, а вот насчет двух АКФ, как это реализовать, мне кажется вот моя ошибка, но я не могу это переписать:
    
    В Matlab во внутреннем цикле (который будет переводом на язык Matlab суммы Matcad) верхний предел в тексте предыдущего сообщения постоянный, а должен быть переменным?
    
    Присваивать надо не RR(i), а RR(k)?
    
    Не Z(i)*Z(i+1), а Z(i)*Z(i+k)?
    
    RR1 = (1./(N-i))*RR; - ?
    - Ответ был удален
      - Н/Д0.00
        
        30.11.2020 11:46
        
        Смотрите, есть спектр сигнала с шумом, есть спектр идеального сигнала, без шума, я смотрю, где спектр полезного сигнала достаточно сильно отличен от нуля, т.е вот эта "горка", а затем я обнуляю составляющие спектральной функции сигнала с шумом, на остальных составляющих, за исключением тех, где спектр полезного сигнала отличен от нуля, достаточно сильно, затем, я должен подвергнуть преобразованию Фурье уже спектральную функцию сигнала полезного с шумом, но с уже обнулены и составляющими, затем сравнить отношение сигнал-шум (мощности) до фильтрации и после
        
        aBoomest+942.89
        
        30.11.2020 11:50
        
        ниже прикрепил.
aBoomest+942.89

30.11.2020 11:49

Давайте сперва со спектром разберемся.

1. Зануление: что а цифры 40 и 80 ? Ваш спектр (приведенный) не имеет таких больших значений. поставил там пока 0,07 на глаз, чтоб видно было результат.

2. Вы почему-то зануляли идеальный спектр, а не спектр с шумом. Может так и надо? Однако поменял.

PS: До строки 100, дальше пока не смотрел. Посмотрите.
- Н/Д0.00
  
  30.11.2020 12:01
  
  Ну последний график не похож на график, который как в 3.png... И просто он строит в нуле две пряме с множеством значение по оси у. Прикладываю скриншот
- Н/Д0.00
  
  30.11.2020 12:02
  
  Я может тупой (вернее не может), но график не строится, скрин приложил, решил убрать subplot, чтобы тоже на одном графике вывести сам график сигнала полезного с шумом до фильтрации, и после
  - aBoomest+942.89
    
    30.11.2020 12:13
    
    Про какой последний? Посли строки 100 в приложенном файле (main2) - не ходил.
    
    Вот графики. Они верные? Если да, то попозже глянем на корреляции.
    
    Отвлеченно, у вас матлаб какой? А то у меня явно с кодировкой русских букв не то. Вот понять бы это после выкачивания с сайта такое или само по себе?
    - Н/Д0.00
      
      30.11.2020 12:15
      
      Они верные, но 5 не не получается, даже если диапазон настроить такой же, как на у Вас на скриншоте...
      - aBoomest+942.89
        
        30.11.2020 12:21
        
        Как это? Т.е. у Вы запускаетет main2 выдает другое? Что выдает? Картинку покажите?
        Вобще учитывая ф-цию рандом, оно от запуска к запуску будет разное. Но должно быть визуально похоже друг на друга.
    - Н/Д0.00
      
      30.11.2020 12:18
      
      Все работает, это у меня на ПК матлаб завис, после перезапуска графики есть! Извините меня!
      - aBoomest+942.89
        
        30.11.2020 12:21
        
        Не проблема. )
      - Н/Д0.00
        
        3.12.2020 18:37
        
        Решил переписать ещё раз часть кода с АКФ, но опять, ошибка, т.е переписать цикл из маткад двумя циклами в матлаб, чтобы получить правильную картинку
        
        % строим две АКФ % чистый шум ZZ1 = Z-Y1; kk = 1:ns/2; RR = zeros(1,(ns/2)); for i = 1:(N-kk-1) for k = 1:ns/2 RR(i) = Z(i)*Z(i+k); end RR1 = (1/(N-i))*RR; end RR1 = zeros(1,(ns/2)); for i = 1:(N-kk-1) for k = 1:ns/2 RR1(i) = ZZ1(i)*ZZ1(i+k); end RR11 = (1/(N-i))*RR1; end figure plot (kk,RR1,kk,RR11); xlim([0 500]);
        
        aBoomest+942.89
        
        3.12.2020 21:32
        
        kk = 1:ns/2;
        RR = zeros(1,(ns/2));
        for i = 1:(N-kk-1)
        for k = 1:ns/2
        RR(i) = Z(i)*Z(i+k);
        end
        RR1 = (1/(N-i))*RR;
        endRR(i) = Z(i)*Z(i+k);
        
        Признаться глубоко "залезть" сегодня возможности не было. Вероятно тут что-то не то.
        
        1. kk - это выходит согласно вашему коду - массив. Тогда вот это for i = 1:(N-kk-1) совнршенно не понятно.
        
        2. RR(i) = Z(i)*Z(i+k); - Z состоит из 4096 элементов. При i+k явно в один прекрасный момент получится больше чем 4096.
        
        ЗЫ: это так, с ходу бросилось в глаза.
        
        ЗЫЗЫ: xcorr чем не подходит? Пол минуты и построил АКФ для идеального, зашумленного и фильтрованного сиганлов. Уж не знаю надо оно или нет . . .
        
        Н/Д0.00
        
        4.12.2020 01:36
        
        На картинке в маткад, там по другому график выглядит, т.е в нуле рост а дальше резкое убывание, для сигнала с шумом, это получается не периодичная функция, а на сигнала просто она в виде косинусоиды. Аналогичную картинку можно получить, если строить АКФ через autocorr (если раскомментировать), но там максимальное значение по амплитуде единица, а на деле оно же другое, xcorr я бы применил, но он не так строит АКФ, как надо для временной фильтрации, чтобы понять, есть сигнал или нет его....
        
        Вообще странно, почему autocorr и xcorr дают разный результат, причем autocorr тот, что ближе к истине
        
        aBoomest+942.89
        
        6.12.2020 22:00
        
        АКФ куска синусойды нагуглите как выглядит. Что у вас в маткаде на картинках - я не могу сказать, но это не АКФ куска синуса - это точно.
        
        Так же не ясно что у вас на картинках с АКФ по оси абсцисс. Там цифра 300. Если это 300 точек, то это мало, т.к. АКФ в 2 р длиннее, т.е. исходник 150 точек. А у вас сигнал имеет 4096 точек. 150 - это маленький кусочек от всего сигнала.
        
        Н/Д0.00
        
        7.12.2020 05:46
        
        Надо построить АКФ по дискретным отчётам (k), как раз то, что в маткаде это правильно, так должны ввглядить эти АКФ, просто через autocorr они по значению оси Y меньше, а вид такой же, а вот через xcorr, такое получить не могу
        
        Н/Д0.00
        
        7.12.2020 05:58
        
        Ну может я ошибаюсь, но, там в файле с методичкой в 11.2 на рис.4 показан вид той АКФ, которая должна получиться, чтобы понять можно обнаружить сигнал или нет
        
        aBoomest+942.89
        
        7.12.2020 14:38
        
        А на вашей картинке
        
        RR - это АКФ синуса с шумом?
        RR1 - это АКФ снус с шумом минус синус = шум?
        
        Н/Д0.00
        
        6.12.2020 17:09
        
        Не получается у меня записать сумму из маткада на языке матлаб, как я ни пытался... как же правильно построить через autocorr, я так же не знаю, к тому же в согласованном фильтре по-прежнему на выходе шум, и нет картинки сигнала, т.е. треугольника в случает в прямоугольного импульса и возрастающей прямой в случае затухающей экспоненты
        
        Н/Д0.00
        
        4.12.2020 01:38
        
        Я просто попытался таким образом переписать сумму из маткад на языке матлаб
    - Н/Д0.00
      
      1.12.2020 11:33
      
      Matlab R2020b, у меня тоже русские буквы поехали
aBoomest+942.89

7.12.2020 14:46

Слева и справа одно и тоже только в разных масштабах.
- Н/Д0.00
  
  7.12.2020 15:24
  
  Вот я и думаю, почему по амплитуде такие огромные значения, и выходит она из середины диапазона 0.8, а как построить АКФ относительно отсчётов к не получается с помощью этой функции, а через autocorr можно, но там тоже с амплитудой немного бардак, обе функции в нуле принимают маскималльное значение, равное единице, а оно в нуле должно по логике принимать значение 2π/ω. Но это когда сдвиг по времени, а вот как сделать сдвиг не по времени, а дискретное, выразился коряво, но я имел ввиду вместо τ, использовать k, и масштаб по оси х, это число точек от 0 до k, как я понял autocorr, как раз строит АКФ по дискретным значениям сдвига, а вот xcorr, я не могу так сделать.
  
  Я Вас понимаю, что утомил Вас, но все очень-очень плохо, потому что больше мне консультироваться просто не с кем, за деньги никто не берёт даже.....
  - aBoomest+942.89
    
    7.12.2020 16:58
    
    1. Не математик, поэтому на счет глубокотеоретического значения амплитуды АКФ не могу сказать. В моей практике ни разу не помню чтоб амплитуда (абсолютное значение) играла какую-то роль. Важна как правило форма.
    
    2. АКФ и ВКФ это ф-ции другой, как бы это сказать, . . . не могу подобрать слово. Что-то в духе domain. Если сигнал - это ф-ция t, то АКФиВКФ это ф-ции взаимного сдвига. Зайдите на вики, там помнится была gifка которая хорошо показывает процесс вычисления.
    
    3. Так вот, то как сдвигать и откуда начинать сдвиг - в математике вероятно это имеет значение и называется по-разному, но в практической деятельности это большого значения не имеет. Поэтому друг относительно друга 2 ф-ции можно сдвигать как душе угодно: вперед, назад, . . .
    
    4. Как следствие у вас есть АКФ (синус ромбик), а ось абсцисс на практике как хочется так и определяйте, как вам удобно. Хотите ноль в центре? Пожалуйста, значит просто взаимный сдвиг вы начали с того, что совместили сигналы, а затем начали сдвигать.
    
    5. Ну и вообще, обычно АКФ и ВКФ нормированы, а след-но амплитуда от 0 до 1. Хотя уверен, что нормирование бывает разное. По мне так вариант от 0 до 1 удобный.
    - Н/Д0.00
      
      7.12.2020 17:13
      
      Про 4 пункт я Вас не совсем понял, что надо сделать, чтобы через xcorr, можно было сделать так, чтобы графики АКФ начинались с нулевой отметки.
      
      Возможно, я не так объяснил, я имел ввиду, что xcorr строит во временной области, а не от числа отсчётов сигнала k, когда проведена его дискретизация...
      - aBoomest+942.89
        
        7.12.2020 19:11
        
        Ф-ция xcorr дает вам ординаты, а массив абсцисс (для рисования графика), как надо так и определяйте.
        PS: Да там есть вариант чтоб ещё lags вернул, но зачем?
        PSPS: Погуглил. Народ просто обрезает, чтоб получить ваш рисунок. Примерно то что и предлагается.
        PSPSPS: Никто не мешает написать свою ф-цию :) и в ней сдвиг делайте так как надо, чтобы насиналось с полного совмещения сигналов.
        PSPSPSPS: И также после вызова xcorr в абсциссу никто не межает влепить не время (временное смещение) а сдвиг в номерах отсчетов.
        
        Н/Д0.00
        
        7.12.2020 19:20
        
        Вот про последние два P.S. я и говорю, что не допираю как написать номально свою функцию, и как прикрутить к xcorr, номера отчётов, буду думать значит как это реализовать, за нормировку, кстати, спасибо!
        
        aBoomest+942.89
        
        8.12.2020 07:51
        
        Так в ссылке в предыдущем посте даже пример есть?! Или я не верно истолковал?
        
        Н/Д0.00
        
        8.12.2020 09:08
        
        Не то, что бы неверно, но это ладно, я все же смог построить требуемую АКФ с вашей помощью!
    - Н/Д0.00
      
      8.12.2020 11:39
      Вот еще не понимаю, вроде записал как надо, пытасюь написать согласованный фильтр, и вот в маткаде картинка правильная, а в матлаб она частично правильная, т.е. при проходе сигнала через фильтр, график такой, что нельзя сказать, обнаруживается сигнал или нет. В случае эксопненциального импульса и прямиугольного импульса, сигнал обнаружен, когда наблюдается картинка похожая на треугольники или "горб", при свертке импульсной переходной характеристики (которая является отражением входного сигнала и смещением в нулевой момент), должны получаться либо треугольник либо горбик, и вот в матлабе никак не могу записать, хотя свертку записл через сумму, циклами, и по идее все должно работать, но я не знаю, что не так...
      
      Там в pdf файле маткад код, чтобы картинками не закидывать...
      
      Код:
      
      % экспоненциально затухающий импульс clear clc Tmax=0.35; T = 0.2; N=4096; dt=Tmax/N; tn=0:dt:Tmax; A=2; % Сигнал % sig=A*(tn>= 0.5)-(A*(tn >= 1)); %+(-A*(tn >= 1)); sig=A*(tn>=T/2).*(exp(-(50*(tn-(T/2))))-(exp(-(250*(tn-(T/2)))).*(tn>=T/2))); % sig_res = sig(tn >=T/2).*(0*sig((tn >= T) & (tn <= Tmax))); % Сигнал + помеха x=sig+1*randn(N+1,1)'; % сигнал с шумом x1 = x-sig; % шум без полезного сигнала L=length(sig); %длина паттерна N=length(x); %длина входного сигнала y=zeros(1,N); y1=zeros(1,N); for i=1:N for p=1:L if (i-p)>0 y(i)=y(i)+x(i-p)*dt; y1(i)=y1(i)+x1(i-p)*dt; % y(i)=y(i-p)+x(p)*dt; % y1(i)=y(i-p)+x1(p)*dt; % y(p)=y(i-p)+x(p)*dt; % y1(p)=y(i-p)+x1(p)*dt; end end end figure(7) % вывод результата subplot (4,1,1); plot(tn,sig);grid;title('signal') % шаблон сигнала subplot (4,1,2); plot(tn,x);grid;title('signal+noise') % реальный сигнал subplot (4,1,3); plot(tn,y*1000);grid;title('filter output') % выход фильтра при свертке с зашумленным сигналом subplot (4,1,4); plot(tn,y1*1000);grid;title('filter output without noise') % выход фильтра при свертке с чистым шумом % прямоугольный импульс clear clc Tmax=0.35; T = 0.2; N=4096; dt=Tmax/N; tn=0:dt:Tmax; A=2; % Сигнал sig=A*(tn >= T/2 )+(-A*(tn >= T)); % Сигнал + помеха x=sig+1*randn(N+1,1)'; % сигнал с шумом x1 = x-sig; % шум без полезного сигнала L=length(sig); %длина паттерна N=length(x); %длина входного сигнала y=zeros(1,N); y1=zeros(1,N); for i=1:N for p=1:L if (i-p)>0 y(i)=y(i)+x(i-p)*dt; y1(i)=y1(i)+x1(i-p)*dt; % y(i)=y(i-p)+x(p)*dt; % y1(i)=y(i-p)+x1(p)*dt; % y(p)=y(i-p)+x(p)*dt; % y1(p)=y(i-p)+x1(p)*dt; end end end figure(8) % вывод результата subplot (4,1,1); plot(tn,sig);grid;title('signal') % шаблон сигнала subplot (4,1,2); plot(tn,x);grid;title('signal+noise') % реальный сигнал subplot (4,1,3); plot(tn,y*1000);grid;title('filter output') % выход фильтра при свертке с зашумленным сигналом subplot (4,1,4); plot(tn,y1*1000);grid;title('filter output without noise') % выход фильтра при свертке с чистым шумом
      - aBoomest+942.89
        
        8.12.2020 13:50
        
        1. Объясните пожалуйста смысл того что вы делаете, не совсем догоняю.
        2. поясните пожалуста строчки цикла
        
        y(i)=y(i)+x(i-p)*dt; y1(i)=y1(i)+x1(i-p)*dt;
        
        Н/Д0.00
        
        8.12.2020 13:53
        
        первая строка получается свертка импульсной хар-ки с сигналом зашумленным, а вторая это свертака импульсной хар-ки с шумом, где нет полезного сигнала
        
        Н/Д0.00
        
        8.12.2020 13:57
        
        я просто не знаю, как иначе сумму записать, т.е саму свертку
        
        aBoomest+942.89
        
        8.12.2020 14:30
        
        А что такое dt?
        
        Н/Д0.00
        
        8.12.2020 14:34
        
        интервал, с которым идет отметка по времени, шаг
Н/Д0.00

8.12.2020 13:53

!
aBoomest+942.89

8.12.2020 15:14
aBoomest+942.89

8.12.2020 15:14
- Н/Д0.00
  
  8.12.2020 20:22
  
  Последний вопросик, а как мне организовать смещение импульсной хар-ки в начало координат? (по шкале временной), просто ее тоже надо выводить на график
aBoomest+942.89

8.12.2020 15:15

4-й график я не понял. В скрипте одно написано, на картинках другое.

А еще видосики коллеги Marat можно на ютубе посмотерть, не мало там по данной теме.
- Н/Д0.00
  
  8.12.2020 15:24
  
  на 4 графике свертка импульсной хар-ки с шумом, без самого полезного сигнала
- Н/Д0.00
  
  8.12.2020 15:41
  
  Marat? Я может путаю, это не MatlabinRussia
  - aBoomest+942.89
    
    8.12.2020 20:12
    
    Да. https://www.youtube.com/watch?v=cRcSiALBfZI&list=PLmu_y3-DV2_kpP8oX_Uug0IbgH2T4hRPL
    - Н/Д0.00
      
      8.12.2020 20:27
      
      Ну да, там 18 уроков вроде, по 5-7 минут, я видел, но не все