Tutorial sobre cálculo de espectros de frequência

Sabemos que um sinal senoidal pode ser escrito de duas formas disitintas:

A forma exponencial é a expressão em Série de Fourier.

Uma representação gráfica do espectro de potência é apresentado a seguir:

SInal senoidal

O script a seguir simula a amostragem de uma onda de tensão senoidal que pode ser descrita através da seguinte equação:

Adota-se:

• 
• 
• 
• 

Escolhe-se uma frequência de amostragem .

Utiliza-se a função fft() para o cálculo do espectro de frequência do sinal.

Gráficos são apresentados dos pontos amostrados no domínio do tempo, do módulo e da fase dos coeficientes da Série de Fourier.

fa=10; % frequencia de amostragem

T=1/fa; % intervalo de amostragem

tfinal=10; % tempo final de referencia

t=0:T:tfinal-T; % vetor de tempo

K=0; % sinal DC adicional

A1=1; % amplitude do sinal

f1=0.5; % frequencia do sinal

phi1=0; % angulo de atraso/avanco de fase

x = K + A1*sin(2*pi*f1*t+phi1);

figure(1);

subplot(3,1,1);

stem(t,x);

xlabel('tempo (s)');

ylabel('Tensão (V)');

grid on;

nx=length(x);

f=0:(fa/(nx-1)):fa;

Cn=fft(x)/nx;

modCn = abs(Cn);

subplot(3,1,2);

stem(f,modCn);

grid on;

xlabel('Frequencia (Hz)');

ylabel('|Cn|');

subplot(3,1,3);

% tol define o valor minimo que abs(y) deve ter para ser considerado

% caso contrario recebe o valor nulo

tol = 1e-1;

Cn(abs(Cn) < tol) = 0;

% calcula a fase de Cn

AngleCn = angle(Cn);

stem(f,AngleCn/pi); % normalizado por pi

grid on

xlabel('Frequencia (Hz)');

ylabel('<Cn/\pi (Rad)')

Observações:

• A teoria demonstra que quando um sinal é amostrado o seu espectro se repete infinitamente a cada múltiplo de tanto no lado positivo quanto no lado negativo,
• A escala de frequências é definida pelo intervalo e contêm o mesmo número de pontos nx das amostras no domínio do tempo,
• Na prática interessa apenas o espectro no intervalo já que o espectro no intervalo corresponde ao espectro no intervalo ,
• A utilização de uma tolerância tol é necessária porque o algoritmo numérico pode retornar valores ínfimos para Cn mas diferentes de zero para a maioria das frequências onde deveria conter valores nulos.
• Para a estimativa do espectro de frequência é necessário a normalização da FFT pelo número de pontos nx,
• O gráfico de fase é normalizado pelo número π para facilitar sua interpretação,
• Obviamente a fase do espectro do sinal depende de , para um sinal senoidal puro devemos fazer .

Espectro de um único lado (One sided spectrum)

Interessa aqui o espectro no intervalo . Utiliza-se o mesmo algoritmo do exemplo anterior selecionando no comando plot() apenas o subrray adequado de modCn e AngleCn.

fa=10; % frequencia de amostragem

T=1/fa; % intervalo de amostragem

tfinal=10; % tempo final de referencia

t=0:T:tfinal-T; % vetor de tempo

K=0; % sinal DC adicional

A1=1; % amplitude do sinal

f1=0.5; % frequencia do sinal

phi1=0; % angulo de atraso/avanco de fase

x = K + A1*sin(2*pi*f1*t+phi1);

figure(2);

subplot(3,1,1);

stem(t,x);

xlabel('tempo (s)');

ylabel('Tensão (V)');

grid on;

nx=length(x);

f=0:(fa/(nx-1)):fa;

Cn=fft(x)/nx;

modCn = abs(Cn);

subplot(3,1,2);

stem(f(1:nx/2+1),modCn(1:nx/2+1));

grid on;

xlabel('Frequencia (Hz)');

ylabel('|Cn|');

subplot(3,1,3);

% tol define o valor minimo que abs(y) deve ter para ser considerado

% caso contrario recebe o valor nulo

tol = 1e-1;

Cn(abs(Cn) < tol) = 0;

% calcula a fase de Cn

AngleCn = angle(Cn);

stem(f(1:nx/2+1),AngleCn(1:nx/2+1)/pi); % normalizado por pi

grid on

xlabel('Frequencia (Hz)');

ylabel('<Cn/\pi (Rad)')

Espectro dos dois lados com frequências simétricas (Double sided spectrum)

O espectro deve ser estimado no intervalo de frequências

Deve ser utilizado o comando do matlab fftshift() como ilustrado abaixo.

fa=10; % frequencia de amostragem

T=1/fa; % intervalo de amostragem

tfinal=10; % tempo final de referencia

t=0:T:tfinal-T; % vetor de tempo

K=0; % sinal DC adicional

A1=1; % amplitude do sinal

f1=0.5; % frequencia do sinal

phi1=0; % angulo de atraso/avanco de fase

x = K + A1*sin(2*pi*f1*t+phi1);

figure(3);

subplot(3,1,1);

stem(t,x);

xlabel('tempo (s)');

ylabel('Tensão (V)');

grid on;

nx=length(x);

f=-fa/2:(fa/(nx-1)):fa/2;

z=fft(x)/nx;

Cn=fftshift(z);

modCn = abs(Cn);

subplot(3,1,2);

stem(f,modCn);

grid on;

xlabel('Frequencia (Hz)');

ylabel('|Cn|');

subplot(3,1,3);

% tol define o valor minimo que abs(y) deve ter para ser considerado

% caso contrario recebe o valor nulo

tol = 1e-1;

Cn(abs(Cn) < tol) = 0;

% calcula a fase de Cn

AngleCn = angle(Cn);

stem(f,AngleCn/pi); % normalizado por pi

grid on

xlabel('Frequencia (Hz)');

ylabel('<Cn/\pi (Rad)')