Aproximação de uma função por soma de senos e cossenos

Exemplos de sinais periódicos.

Algumas vezes, desejamos aproximar um sinal periódico (possivelmente ruidoso) por uma soma de senos e cossenos. Sabemos que (quase) qualquer sinal periódico, por exemplo uma onda quadrada, pode ser decomposto em uma soma infinita de senos e cossenos cujas frequências são múltiplas da frequência fundamental do sinal periódico. Essa é a famosa Série de Fourier (ver aqui).

Nesta postagem, usando o Scilab e o comando fft, vamos ver como encontrar as principais componentes de seno e cosseno de um sinal periódico arbitrário. A ideia básica é calcular o espectro do sinal periódico com o comando fft e, em seguida, "montar" um sinal com as componentes de maior amplitude. Essa estratégia, se o sinal tiver vários períodos, acaba por filtrar algum ruído existente na função.

Seja o sinal indicado abaixo, quais suas componentes em termos de seno e cosseno?

Na figura acima, temos 10 períodos de um sinal periódico ruidoso, seu período vale 100 unidades de tempo (= 1 segundo). A análise é feita com o comando fft:

Com essas componentes, chegamos ao seguinte sinal: $$s(t) = 1.009\cos(2\pi t) + 0.495\sin(4\pi t) - 0.3013\cos(6\pi t) + 0.046\sin(8\pi t) + 0.0831\cos(8\pi t) + 0.052\cos(10\pi t)$$.
Para comparação, o sinal "exato" é:
$$s(t) = \cos(2\pi t) + 0.5\sin(4\pi t) - 0.3\cos(6\pi t) + + 0.1\cos(8\pi t + 0.5) + 0.05\cos(10\pi t)$$. Graficamente:

Código Scilab:

close; close;
t=0:0.01:(10-0.01);
p2=2*%pi;
s = cos(p2*t) + 0.5*sin(2*p2*t) - 0.3*cos(3*p2*t) +...
 0.1*cos(4*p2*t + 0.5) + 0.05*cos(5*p2*t);
sp = s;
s = s + 0.1*rand(s,'n');
plot(s); xgrid;

f = 1:max(size(s));
f = f - 1;
f = f/max(f); f = f*100;
sf = fft(s);
sfr = real(sf);
sfc = imag(sf);
// subplot(2,1,2);
figure;
plot(f(1:200), sfr(1:200),'m');
plot(f(1:200), sfc(1:200),'r');
legend('Componentes em fase','Componentes em quadratura');
//pega picos:
tt = 1:max(size(s));
tt=tt-1; tt=tt*0.01;
p = zeros(1,5);
fp = p;
ss = 0*tt;
tam2 = round(max(size(sfr))/2);
sfr = sfr(1:tam2); sfr(1) = 0;
sfc = sfc(1:tam2); sfc(1) = 0;
m = mean(s);
picoc = 0
picos = 0;
kf = 500; //fator de escala - 10 períodos e dt = 0,01;
for k=1:6
    [a,b] = max(abs(sfr));
    fq = round(f(b));
    if abs(sfr(b))>picoc then picoc = abs(sfr(b)); end;
    if abs(sfr(b))>(picoc/30) then
        ss = ss + sfr(b)*cos(fq*p2*tt)/kf;
        disp([fq,sfr(b)/kf,0]);
    end;
    sfr(b) = 0;
    [a,b] = max(abs(sfc));
    fq = round(f(b));
    if abs(sfc(b))>picos then picos = abs(sfc(b)); end;  
    if abs(sfc(b))>(picoc/30) then
        ss = ss - sfc(b)*sin(fq*p2*tt)/kf;
        disp([fq,sfc(b)/kf,1]);
    end;
    sfc(b) = 0;
end
ss = m + ss;
figure; plot(tt,sp,'m',tt,ss,'b',tt,s,'r');
legend('Sinal puro','Sinal calculado','Sinal com ruído');