Densidade espectral de potência usando o Scilab

A densidade espectral de potência (DEP ou power spectral density - PSD) é uma função positiva que indica como está distribuída a energia de um sinal em função da frequência. Se o sinal for periódico, a DEP será composta por impulsos concentrados nas frequências que formam esse sinal. Nos casos mais gerais, o sinal pode ser um processo aleatório e o cálculo da sua densidade espectral de potência é bem mais complicada.

O cálculo da DEP é importante, por exemplo, quando se trabalha com códigos de linha ou em problema de filtragem. Para um sinal (código de linha) formado por uma sequência de pulsos, DEP pode ser calculada como segue.

Sistema que tem por saída os pulsos y(t) (código de linha).

A densidade espectral de potência de um código de linha pode ser calculada por:
$$ S_y(f) = |P(f)|^2 S_x(f) $$
onde $P(f)$ é transformada de Fourier do pulso $p(t)$ e $S_x(f)$ é a densidade espectral de potência do trem de impulsos que gera o sinal $y(t)$. O sinal $x(t)$ depende da sequência de bits e da codificação usada e sua função autocorrelação pode ser calculada por:
$$ R_n = \lim_{N \rightarrow \infty} \frac{1}{N}\sum_{k}a_k a_{k+n} = E(a_k a_{k+n}) $$
e
$$ S_x(f) = \frac{1}{T_b}\sum_{n = -\infty}^{n = \infty} R_n e^{-jn2\pi f T_b} $$
Logo
$$ S_y(f) = \frac{|P(f)|^2}{T_b} \left[\sum_{n = -\infty}^{n = \infty} R_n e^{-jn2\pi f T_b} \right] $$
ou ainda
$$ S_y(f) = \frac{|P(f)|^2}{T_b} \left[R_0 + 2\sum_{n = 1}^{n = \infty} R_n cos(2\pi n f T_b) \right] $$
onde
$$ R_0 = \lim_{N \rightarrow \infty} \frac{1}{N}\sum_{k}a_k^2 = E(a_k^2). $$
Por exemplo, para uma sinalização polar RZ (retorno ao zero), é possível mostrar que
$$ S_y(f) = \frac{T_b}{4} sinc^2\left( \frac{\pi f T_b}{2} \right) $$

Toda essa matemática pode ser simulada usando-se o Scilab. Para o exemplo acima, obtemos o seguinte gráfico com o código abaixo.

Pulsos aleatórios (uma realização), DEP obtida pela simulação e curva teórica.

Código:

Np = 500; // número de pulsos
mmx = 199;  // número de curvas para obter a média

// pulso p(t) 1:
p = [ones(1,50), zeros(1,50)];
pf = fft(p);
pf2 = abs(pf);
subplot(2,2,1); plot(pf2);

// sinal aleatório
ak = sign(rand(1,1,'n'));
sinal = ak*p;
for k=1:(Np-1)
    ak = sign(rand(1,1,'n'));
    sinal = [sinal, ak*p];
end
subplot(2,2,2); plot(sinal*1.01);
fss = fft(sinal);
fss2 = abs(fss).*abs(fss);

fsk = fss2;

for mm=1:mmx
sinal = sign(rand(1,1,'n'))*p;
for k=1:(Np-1)
    ak = sign(rand(1,1,'n'));
    sinal = [sinal, ak*p];
end
fss = fft(sinal);
fss2 = abs(fss).*abs(fss);
fsk = fsk + fss2;
end;

fsk = fsk/(mm+1);
fsk = fsk/(Np*50*50*4);
subplot(2,2,3); plot(fsk);

// Sy teorico:
f = 1:1000; 
f = f/100;
sy1 = 0.25*(sin(%pi*f)./(%pi*f)).^2;
plot(f*1000,sy1,'m');