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

[Alm99]

Добрый вечер, уважаемые форумчане! Прошу помощи, необходимо в системе Matlab повторить все аналогичные действия, перечисленные в методичке пункты 2-5.

Проблема в том, что не могу построить обратное преобразование Фурье и вывести его, а также построить корреляционные (АКФ) функции сигналов. Также при записи в согласованной фильтрации, при ее выполнении для экспоненциально затухающего импульса и прямоугольного импульса, начинающихся в середине интервала, на выходе фильтра получается не то что надо (для экспоненциально затухающего импульса на выходе должна быть линия, возрастающая после начала импульса в середине, а для прямоугольного импульса должен получиться график в виде треугольника, максимум которого достигается в точке начала следования прямоугольного импульса).
Я повторил аналогичные действия в Mathcad, чтобы визуализировать то, что надо сделать в матлабе. На всякий случай скину методичку и сами скриншоты из маткада. Не получается повторить все те же действия, что и в маткаде, но с помощью матлаб... Не могу сделать, даже за деньги заказ не берут, а осталось всего-то, но вот не могу исправить ошибки, которые выдает программа! Буду благодарен за оказанную помощь!







Методические указания : https://drive.google.com/file/d/16MEhu72rEczrhVsHcTUV5OXZCMST1EN-/view?usp=sharing

Вот код в матлаб:

Используется синтаксис Matlab M

% % узкополосная фильтрация
% 
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 200])

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


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

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


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

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


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

[Alm99]

Немного подправил, но все равно, график не строится, ни для АКФ, ни для обратного преобразования Фурье, ни для согласованного фильтра

Используется синтаксис Matlab M

% % узкополосная фильтрация
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 100])

% обнуляем те элменты к спектральной функции...
%  к, где спектр полезного сигнала  не  сильно отличается от нуля, и
%  совершаем обратное преобразование Фурье
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,(N-ns-1));
for i = 1:N-ns-1
    RR(i) = Z(i)*Z(i+1);
end
RR = (1./(N-i)).*RR;


RR1 = zeros(1,(N-ns-1));
for i = 1:N-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') % выход фильтра без учета шума

-- 24.11.2020, 23:24 --

не форматируется код почему-то

Pphantom:

Потому что не надо вставлять тэги подсветки синтаксиса внутрь тэга кода, они независимы. Поправил.

[GAA]

В тексте предыдущего сообщения есть фрагмент

Используется синтаксис Matlab M

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

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

N = 4096, ns = 4096. На строчке RR = (1./(N-i)).*RR; завершается выполнение с сообщением об ошибке:

??? Error using ==> .*
Matrix dimensions must agree.

Error in ==> D:\inet\dxdy\MatLab\Alm99.m
On line 79  ==> RR = (1./(N-i)).*RR;

Т.е. Matlab 6.5 указывает на место с проблемой.
Upd. Matlab 2013 сообщает аналогично

> In Alm99 at 69 
Error using  .* 
Matrix dimensions must agree.

Error in Alm99 (line 79)
RR = (1./(N-i)).*RR;

[Alm99]

GAA


Ну теперь все равно, я не понимаю, почему он не хочет строить две АКФ

Используется синтаксис Matlab M

% строим две АКФ
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);
 

[GAA]

Matlab 2013 сообщает

Index exceeds matrix dimensions.

Error in Alm99 (line 79)
    RR(i) = Z(i)*Z(i+1);

И это ожидаемо: вылезаем за границу (см. число элементов и максимальный индекс элемента к которому идёт обращение в цикле)

[Alm99]

GAA


Но я все равно не могу построить АКФ

Используется синтаксис Matlab M

ZZ1 = Z-Y1;

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


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

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

-- 25.11.2020, 17:50 --

GAA
И самое главное, я не могу понять, почему свертка работает не так, как надо... Почему-то результат получается не тот, что должен быть. неделю не могу написать скрипт, хотя ошибка наверное тривиальная

[GAA]

За блоком RR идёт блок RR1

Используется синтаксис Matlab M

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

Matlab сообщает об ошибке

Error using  .* 
Matrix dimensions must agree.

Error in Alm99 (line 89)
RR1 = (1./(N-i)).*RR1;

В следующий раз, пожалуйста, читайте предупреждения и сообщения Matlab об ошибках.
И если сами не можете исправить, то сообщайте о них [ошибках и предупреждениях] в сообщении на форуме. А если исправили, то указывайте всё что исправили.

[Alm99]

GAA

Простите меня, уважаемый GAA! Но сейчас нет такой ошибки, но график все равно не строится, я не могу поправить ошибку....

Используется синтаксис Matlab M

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

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


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

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

Т.е не получается повторить то же, как и в маткаде... Сейчас нет ошибок, как и нет графика, хотя там должен быть сдвиг, в АКФ и я возможно задал его неправильно...

[GAA]

У Вас в Mathcad k — переменная диапазон. При помощи этой переменной вы создаёте матрицы RR_k и RR1_k.

В тексте Matlab ns у Вас — скаляр: ns = (length(tt));

Естественно, циклы вычисления RR и RR1 нужно переделать: верхний предел сумм в тексте Mathcad не постоянная величина, а зависит от номера элемента матрицы (RR или RR1). До этого места там ещё места есть странные.

[Alm99]

GAA


Если Вы про то, как записано в маткаде, то это правильно, все как по пунктам в методичке, ну а в Matlab так сумму записать не получается, можно так, но это не две АКФ:

Используется синтаксис Matlab M

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 (5)
plot (kk,RR1,kk,RR11);
xlim([0 500]);

[GAA]

Я про то, как записано на Matlab.
В Mathcad двойной цикл (внешний по k, а внутренний по i). Можно дословно перевести на Matlab: внешний цикл в Matlab сделать по k, а внутренний по i ("во внешнем цикле менять верхний предел внутреннего".)

[Alm99]

GAA
Я конечно не понял как это сделать, но попытаюсь...

[Alm99]

GAA


Все равно не выходит, решил еще попробовать через xcorr(), тоже не то строит, что нужно...

Используется синтаксис Matlab M

% строим две АКФ
ZZ1 = Z-Y1; % чистый шум

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


RR1 = zeros(1,(ns/2));
for k = 1:ns/2
    for i = 1:ns/2
    RR1(i) = ZZ1(i)*ZZ1(i+1);
    end
    k = k-i;
end
RR11 = (1./(N-i))*RR1;
plot (kk,RR1,kk,RR11);
xlim([0 500]);
ylim([-2 8]);

% Corr_auto1 = xcorr(Z);
% Corr_auto2 = xcorr(ZZ1);
% figure (6)
% plot (1:length(Corr_auto1),Corr_auto1,1:length(Corr_auto2),Corr_auto2);
% xlim([0 500]);
 

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

[GAA]

Вычисление RR в Mathcad из начального сообщения темы:

Вложение:

Alm_Mathcad_tf.PNG

1. В Matlab во внутреннем цикле (который будет переводом на язык Matlab суммы Matcad) верхний предел в тексте предыдущего сообщения постоянный, а должен быть переменным.
2. После внутреннего цикла изменять переменную k не надо, т.е. k=k-1 в тексте надо убрать.
3. Присваивать надо не RR(i), а RR(k).
4. Не Z(i)*Z(i+1), а Z(i)*Z(i+k).
5. RR1 = (1./(N-i))*RR; - очевидно неправильно.

Кроме этого, в Mathcad переменная k имеет другой диапазон значений (от 0 до N/2).
Перевод с Mathcad на Matlab прост. Поэтому если в тексте Mathcad нет ошибок, то написать текст легко.

[Alm99]

GAA в сообщении #1494227 писал(а):

1. В Matlab во внутреннем цикле (который будет переводом на язык Matlab суммы Matcad) верхний предел в тексте предыдущего сообщения постоянный, а должен быть переменным.

Т.е получается верхний предел есть сам массив еще? Или он будет как ns-k?

Получилось построить вот так, но значения по оси Y, меньше...

Используется синтаксис Matlab M

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]);
 

А если делать предел верхний не постоянным, то так?

Используется синтаксис Matlab M

%строим две АКФ
ZZ1 = Z-Y1; % чистый шум

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

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

GAA:

Или Вы издеваетесь над доброжелательностью форума, или Вы не пытаетесь изучить Matlab, но почему-то переводите на Matlab текст рабочего листа Mathcad. Если издеваетесь, то перестаньте. А если ещё не начали изучать Matlab, то откройте рекомендованный Вам учебник. Если никакой учебник не рекомендовали, то можно посмотреть

GAA в сообщении #1447078 писал(а):

Ануфриев И. Е., Смирнов А. Б., Смирнова Е. Н. MATLAB 7.
Н.Н.Мартынов, А.П.Иванов MATLAB 5.X. Вычисления, визуализация, программирование.

Это не форум для выполнения учебных лабораторных работ за студентов или преподавателей (и студенты, и преподаватели должны учебные работы в основном делать сами).
Две недели на изучение Matlab.