2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему
 
 Решить задачу обнаружения сигнала в системе Matlab
Сообщение24.11.2020, 19:08 


17/03/20
183
Добрый вечер, уважаемые форумчане! Прошу помощи, необходимо в системе 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') % выход фильтр
 

 Профиль  
                  
 
 Re: Решить задачу обнаружения сигнала в системе Matlab
Сообщение24.11.2020, 23:24 


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

код: [ скачать ] [ спрятать ]
Используется синтаксис 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 --

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

 Профиль  
                  
 
 Re: Решить задачу обнаружения сигнала в системе Matlab
Сообщение25.11.2020, 16:59 
Заслуженный участник


12/07/07
4438
В тексте предыдущего сообщения есть фрагмент
Используется синтаксис 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;

 Профиль  
                  
 
 Re: Решить задачу обнаружения сигнала в системе Matlab
Сообщение25.11.2020, 17:26 


17/03/20
183
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);
 

 Профиль  
                  
 
 Re: Решить задачу обнаружения сигнала в системе Matlab
Сообщение25.11.2020, 17:35 
Заслуженный участник


12/07/07
4438
Matlab 2013 сообщает
Код:
Index exceeds matrix dimensions.

Error in Alm99 (line 79)
    RR(i) = Z(i)*Z(i+1);
И это ожидаемо: вылезаем за границу (см. число элементов и максимальный индекс элемента к которому идёт обращение в цикле)

 Профиль  
                  
 
 Re: Решить задачу обнаружения сигнала в системе Matlab
Сообщение25.11.2020, 17:46 


17/03/20
183
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
И самое главное, я не могу понять, почему свертка работает не так, как надо... Почему-то результат получается не тот, что должен быть. неделю не могу написать скрипт, хотя ошибка наверное тривиальная

 Профиль  
                  
 
 Re: Решить задачу обнаружения сигнала в системе Matlab
Сообщение25.11.2020, 18:05 
Заслуженный участник


12/07/07
4438
За блоком 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 об ошибках.
И если сами не можете исправить, то сообщайте о них [ошибках и предупреждениях] в сообщении на форуме. А если исправили, то указывайте всё что исправили.

 Профиль  
                  
 
 Re: Решить задачу обнаружения сигнала в системе Matlab
Сообщение25.11.2020, 18:15 


17/03/20
183
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);


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

 Профиль  
                  
 
 Re: Решить задачу обнаружения сигнала в системе Matlab
Сообщение25.11.2020, 21:53 
Заслуженный участник


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

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

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

 Профиль  
                  
 
 Re: Решить задачу обнаружения сигнала в системе Matlab
Сообщение25.11.2020, 22:16 


17/03/20
183
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]);

 Профиль  
                  
 
 Re: Решить задачу обнаружения сигнала в системе Matlab
Сообщение25.11.2020, 22:39 
Заслуженный участник


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

 Профиль  
                  
 
 Re: Решить задачу обнаружения сигнала в системе Matlab
Сообщение26.11.2020, 00:03 


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

 Профиль  
                  
 
 Re: Решить задачу обнаружения сигнала в системе Matlab
Сообщение26.11.2020, 20:07 


17/03/20
183
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]);
 


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

 Профиль  
                  
 
 Re: Решить задачу обнаружения сигнала в системе Matlab
Сообщение26.11.2020, 22:36 
Заслуженный участник


12/07/07
4438
Вычисление RR в Mathcad из начального сообщения темы:
Вложение:
Alm_Mathcad_tf.PNG
Alm_Mathcad_tf.PNG [ 17.22 Кб | Просмотров: 1960 ]

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 нет ошибок, то написать текст легко.

 Профиль  
                  
 
 Re: Решить задачу обнаружения сигнала в системе Matlab
Сообщение26.11.2020, 23:21 


17/03/20
183
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.

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 15 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group