quarta-feira, 9 de dezembro de 2015

Scilab: procurando por picos em um sinal



Em algumas situações estamos interessados em encontrar os valores máximos (ou mínimos) em um sinal qualquer. Nesta postagem mostramos duas formas de conseguirmos isso. Na primeira, é feita uma busca exaustiva por picos positivos e negativos depois de se estabelecer um limiar de decisão*. Na segunda, é usado simplesmente o comando "gsort" do Scilab. Em seguida, é mostrado o sinal "janelado" em torno dos picos significativos. Segue o código com alguns poucos cometários.

* é muito difícil de se estabeler esse limiar sem conhecimento prévio de algumas características do sinal a ser analisado.


Código

///// Encontrando picos em um sinal
///// Prof. Francisco José
///// 09/12/2015
///// Versão 01.

close; close; close;
clear; clc;

pi = %pi;

// Gerando um sinal com picos:
fa = 12; dt = 1/fa;
t = 0:dt:5; w1 = 2*pi*1.5; w2 = 2*pi*1.8;
s = 2*sin(w1*t) + 2*cos(w2*t);
s = s + 0.5*rand(1,max(size(t)),'n');
p = [0:dt:2.5, 2.5:-dt:dt];
s = s.*p; s = s + 2*sin(4*w1*t);
s = [s, 0.1*s, -s, 0*s, 0.1*s, s, -0.1*s];
h = [0.2 0.8 0.9 0.6 0.3 0.1];
s = filter(h,1,s);
s = s + 0.5*rand(1,max(size(s)),'n');
s = s + filter([0.7 0.2 0.1],1,s);

// Sinal com picos e ruido gerado:
plot(s);

// Procurando picos positivos e negativos:
sa = abs(s);
smedio = mean(sa);
smaxi = max(sa);
razao = smaxi/smedio;
disp(razao);

if razao < 2 then  disp('Não existem picos significativos');
end
limiar = 0.55*(smaxi + smedio);


vpk = [];  // vetor com os locais dos picos máximos
k = 1;
while (k < max(size(s)))
    // picos positivos:
    if s(k) > limiar then
        while s(k+1)>s(k)
            k = k + 1;
        end
        vpk = [vpk, k];
    end
    // picos negativos:
    if s(k) < (-limiar) then
        while (s(k+1) < s(k))
            k = k + 1;
        end
        vpk = [vpk, k];
    end
    k = k + 1;
end

figure; plot(s);
plot(vpk,s(vpk),'ro',[1 k],[limiar, limiar],[1 k],[-limiar, -limiar]);
title('Algoritmo de busca de picos');

////////// Usando o gsort:
[skk,kk] = gsort(sa);
vpk = kk(1:12);  // os 12 maiores picos
figure; plot(s); plot(vpk,s(vpk),'ro'); xgrid;
title('Usando o gsort: 12 maiores picos');

//// Isolando os picos:
mk = 0*s;
for k=1:10
    mk(vpk(k)-5:vpk(k)+5) = 1;
end
mk = convol(mk,[0.1 0.5 0.3 0.1]); // suavizando a máscara
mk = mk(1:max(size(s)));
s_iso = s.*mk;
figure; plot(s_iso); title('Regiões com picos isolados');
plot(mk,'r');

Figuras:


Nenhum comentário:

Postar um comentário