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:
fftseno.png
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:
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)')