• Регистрация
© 2020г. Веб Инновации
  • +0
Н/Д
Н/Д 0.00
н/д

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

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') % выход фильтра без учета шума

Теги

      25.11.2020
      • 55
      • 1196

      Комментарии

      • aBoomest
        aBoomest+942.89
        26.11.2020 05:40

        1. А картинки в маткаде априори правильные?
        2. До БПФ поправил. Дальше смотрите, правильнные картинки или как? Если нет, то желательно конкретные вопросы.

      • aBoomest
        aBoomest+942.89
        26.11.2020 09:07

        Не, все должно работать. Скольки процентные (отн. макс.) значения вы обнуляете?

      • aBoomest
        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 - это у вас почему то модуль спектра, а не сам спектр. Вы модуль берите потом когда рисуете график, либо промежуточную переменную дополнительную надо. А сам спектр - это совсем не модуль.

        • Ой, но так тогда у меня не получается переписать, чтобы сначала был спектр, затем нанести его модуль на график, а потом уже вычленить незначащие частоты и построить результат ОБПФ вместе с исходным зашумленным сигналом, т.е. как на скринах в маткаде

          • aBoomest
            aBoomest+942.89
            30.11.2020 10:35

            Немного не так.

            1. У вас д.б. спектр - как есть. Спектр - есть спектр.

            2. По амплитуде спектра вы должны определить элементы массива подлежащие "занулению".

            3. И занулить их надо именно в спектре, а не в модуле, иначе от чего делать обратную процедуру.
            А отрисовывать конечно надо либо амплитуду и фазу, либо re и im, это уж как удобно.

            • Попытаюсь конечно переписать, хотя  3 недели уже никак не могу сделать, а вот насчет двух АКФ, как это реализовать, мне кажется вот моя ошибка, но я не могу это переписать:

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


              Присваивать надо не RR(i), а RR(k)?

              Не Z(i)*Z(i+1), а Z(i)*Z(i+k)?

              RR1 = (1./(N-i))*RR; - ?

              • Ответ был удален

                • Смотрите, есть спектр сигнала с шумом, есть спектр идеального сигнала, без шума, я смотрю, где спектр полезного сигнала достаточно сильно отличен от нуля, т.е вот эта "горка", а затем я обнуляю составляющие спектральной функции сигнала с шумом, на остальных составляющих, за исключением тех, где спектр полезного сигнала отличен от нуля, достаточно сильно, затем, я должен подвергнуть преобразованию Фурье уже спектральную функцию сигнала полезного с шумом, но с уже обнулены и составляющими, затем сравнить отношение сигнал-шум (мощности)  до фильтрации и после

                  • aBoomest
                    aBoomest+942.89
                    30.11.2020 11:50

                    ниже прикрепил.

          • aBoomest
            aBoomest+942.89
            30.11.2020 11:49

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

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

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

            PS: До строки 100, дальше пока не смотрел. Посмотрите.

            • Я может тупой (вернее не может), но график не строится, скрин приложил, решил убрать subplot, чтобы тоже на одном графике вывести сам график сигнала полезного с шумом до фильтрации, и после

              • aBoomest
                aBoomest+942.89
                30.11.2020 12:13

                Про какой последний? Посли строки 100 в приложенном файле (main2) - не ходил.

                Вот графики. Они верные? Если да, то попозже глянем на корреляции.

                Отвлеченно, у вас матлаб какой? А то у меня явно с кодировкой русских букв не то. Вот понять бы это после выкачивания с сайта такое или само по себе?

                • Они верные, но 5 не не получается, даже если диапазон настроить такой же, как на у Вас на скриншоте...

                  • aBoomest
                    aBoomest+942.89
                    30.11.2020 12:21

                    Как это? Т.е. у Вы запускаетет main2 выдает другое? Что выдает? Картинку покажите?
                    Вобще учитывая ф-цию рандом, оно от запуска к запуску будет разное. Но должно быть визуально похоже друг на друга.

                  • Все работает, это у меня на ПК матлаб завис, после перезапуска графики есть! Извините меня!

                    • aBoomest
                      aBoomest+942.89
                      30.11.2020 12:21

                      Не проблема. )

                       

                      • Решил переписать ещё раз часть кода с АКФ, но опять, ошибка, т.е переписать цикл из маткад двумя циклами в матлаб, чтобы получить правильную картинку

                         

                        % строим две АКФ

% чистый шум 
 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
                          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 чем не подходит? Пол минуты и построил АКФ для идеального, зашумленного и фильтрованного сиганлов. Уж не знаю надо оно или нет . . .

                          • На картинке в маткад, там по другому график выглядит, т.е в нуле рост а дальше резкое убывание, для сигнала с шумом, это получается не периодичная функция, а на сигнала просто она в виде косинусоиды. Аналогичную картинку можно получить, если строить АКФ через autocorr (если раскомментировать), но там максимальное значение по амплитуде единица, а на деле оно же другое, xcorr я бы применил, но он не так строит АКФ, как надо для временной фильтрации, чтобы понять, есть сигнал или нет его....

                             

                            Вообще странно, почему autocorr и xcorr дают разный результат, причем autocorr тот, что ближе к истине

                             

                            • aBoomest
                              aBoomest+942.89
                              6.12.2020 22:00

                              АКФ куска синусойды нагуглите как выглядит. Что у вас в маткаде на картинках - я не могу сказать, но  это не АКФ куска синуса - это точно.

                              Так же не ясно что у вас на картинках с АКФ по оси абсцисс. Там цифра 300. Если это 300 точек, то это мало, т.к. АКФ в 2 р длиннее, т.е. исходник 150 точек. А у вас сигнал имеет 4096 точек. 150 - это маленький кусочек от всего сигнала.

                            • Ну может я ошибаюсь, но, там в файле с методичкой в 11.2 на рис.4 показан вид той АКФ, которая должна получиться, чтобы понять можно обнаружить сигнал или нет

                              • aBoomest
                                aBoomest+942.89
                                7.12.2020 14:38

                                А на вашей картинке

                                RR - это АКФ синуса с шумом?
                                RR1 - это АКФ снус с шумом минус синус = шум?

                • aBoomest
                  aBoomest+942.89
                  7.12.2020 14:46

                  Слева и справа одно и тоже только в разных масштабах.

                  • Вот я и думаю, почему по амплитуде такие огромные значения, и выходит она из середины диапазона 0.8, а как построить АКФ относительно отсчётов к не получается с помощью этой функции, а через autocorr можно, но там тоже с амплитудой немного бардак, обе функции в нуле принимают маскималльное значение, равное единице, а оно в нуле должно по логике принимать значение 2π/ω. Но это когда сдвиг по времени, а вот как сделать сдвиг не по времени, а дискретное, выразился коряво, но я имел ввиду вместо τ,  использовать k,  и масштаб по оси х, это число точек от 0 до k,  как я понял autocorr, как раз строит АКФ по дискретным значениям сдвига, а вот xcorr, я не могу так сделать.

                     

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

                     

                     

                     

                    • aBoomest
                      aBoomest+942.89
                      7.12.2020 16:58

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

                      2. АКФ и ВКФ это ф-ции другой, как бы это сказать,  . . . не могу подобрать слово. Что-то в духе domain. Если сигнал - это ф-ция t, то АКФиВКФ это ф-ции взаимного сдвига. Зайдите на вики, там помнится была gifка которая хорошо показывает процесс вычисления.

                      3. Так вот, то как сдвигать и откуда начинать сдвиг - в математике вероятно это имеет значение и называется по-разному, но в практической деятельности это большого значения не имеет. Поэтому друг относительно друга 2 ф-ции можно сдвигать как душе угодно: вперед, назад, . . .

                      4. Как следствие у вас есть АКФ (синус ромбик), а ось абсцисс на практике как хочется так и определяйте, как вам удобно. Хотите ноль в центре? Пожалуйста, значит просто взаимный сдвиг вы начали с того, что совместили сигналы, а затем начали сдвигать.

                      5. Ну и вообще, обычно АКФ и ВКФ нормированы, а след-но амплитуда от 0 до 1. Хотя уверен, что нормирование бывает разное. По мне так вариант от 0 до 1 удобный.

                      • Про 4 пункт я Вас не совсем понял, что надо сделать, чтобы через xcorr, можно было сделать так, чтобы графики АКФ начинались с нулевой отметки.

                         

                         Возможно, я не так объяснил, я имел ввиду, что xcorr строит во временной области, а не от числа отсчётов сигнала k, когда проведена его дискретизация...

                        • aBoomest
                          aBoomest+942.89
                          7.12.2020 19:11

                          Ф-ция xcorr дает вам ординаты, а массив абсцисс (для рисования графика), как надо так и определяйте.
                          PS: Да там есть вариант чтоб ещё lags вернул, но зачем?
                          PSPS: Погуглил. Народ просто обрезает, чтоб получить ваш рисунок. Примерно то что и предлагается.
                          PSPSPS: Никто не мешает написать свою ф-цию :) и в ней сдвиг делайте так как надо, чтобы насиналось с полного совмещения сигналов.
                          PSPSPSPS: И также после вызова xcorr в абсциссу никто не межает влепить не время (временное смещение) а сдвиг в номерах отсчетов.

                          • Вот про последние два P.S. я и говорю, что не допираю как написать номально свою функцию, и как прикрутить к xcorr, номера отчётов, буду думать значит как это реализовать, за нормировку, кстати, спасибо!

                            • aBoomest
                              aBoomest+942.89
                              8.12.2020 07:51

                              Так в ссылке в предыдущем посте даже пример есть?! Или я не верно истолковал?

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

                        Там в 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
                          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;

                          • я просто не знаю, как иначе сумму записать, т.е саму свертку

                            • aBoomest
                              aBoomest+942.89
                              8.12.2020 14:30

                              А что такое dt?

                • aBoomest
                  aBoomest+942.89
                  8.12.2020 15:14
                  • aBoomest
                    aBoomest+942.89
                    8.12.2020 15:14
                  • aBoomest
                    aBoomest+942.89
                    8.12.2020 15:15

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

                    А еще видосики коллеги Marat можно на ютубе посмотерть, не мало там по данной теме.

                    • Marat? Я может путаю, это не MatlabinRussia

                      • aBoomest
                        aBoomest+942.89
                        8.12.2020 20:12

                        Да. https://www.youtube.com/watch?v=cRcSiALBfZI&list=PLmu_y3-DV2_kpP8oX_Uug0IbgH2T4hRPL

                  Популярные посты

                  Темы

                  Популярные теги

                  MATLAB (305) Simulink (201) конкурс_hub.exponenta_2019 (54) конкурс_hub_2021 (46) Student_Hub_2021 (45) Электропривод (36) StudentHub2019 (35) системы управления (35) VIIВСЕРОССИЙСКАЯКОНФЕРЕНЦИЯ (30) доклад (28) моделирование (27) математическое моделирование (26) MATLAB-программирование (25) электротехника (25) МОП (25)
                  В скобках указано количество постов по тегу
                  Все теги