Universidade Federal de Minas Gerais

Departamento de Engenharia Eletrônica

ELT129 – Oficina de Modelagem e Simulação

Professor: Leonardo Mozelli – lamoz@[Link]

Tutorial 9 – Análise Frequencial

1 Introdução
Este tutorial tem o intuito de familiarizar o leitor com os Diagramas de Bode e a ferramenta
conhecida como FFT (Fast Fourier Transform).

2 Fast Fourier Transform – FFT

FFT é o nome dado a algoritmos rápidos que calculam a DFT (Discrete Fourier Transform).
Na DFT, o sinal é discreto nos domı́nios do tempo e da frequência e, portanto, pode ser
representado em um computador digital. Note a diferença com relação à DTFT (Discrete
Time Fourier Transform), que transforma um sinal discreto no tempo em um sinal contı́nuo
na frequência.
Uma das propriedades da DTFT afirma que a um sinal discreto no tempo corresponde
um sinal periódico na frequência. Analogamente, um sinal discreto no domı́nio da frequência
implica em um sinal periódico no tempo. Dessa forma, pode-se entender a DFT como sendo
a DTFT de um sinal discreto e finito no tempo que foi estendido para um sinal discreto e
periódico na frequência.
Uma DFT mapeia um sinal discreto no tempo x[n], n = 0, 1, 2, . . . , N − 1 para um sinal
discreto na frequência X[ω]. À sequência temporal x[n],
x[0], x[1 · Ts ], x[2 · Ts ], . . . , x[(N − 1) · Ts ],
corresponde uma sequência frequencial X[ω]:
2πFs 2πFs 2πFs
X[0], X[1 · ], X[2 · ], . . . , X[(N − 1) · ],
N N N
sendo Ts o tempo de amostragem, Fs = 1/Ts a frequência de amostragem e N o número de
pontos disponı́veis. Aqui há duas observações importantes:
• o número de pontos da DFT é o mesmo do sinal original;
• a resolução no domı́nio da frequência (2πFs /N ) será tanto melhor quanto maior o número
N de pontos amostrados para uma taxa de amostragem fixa.
Ainda é importante notar que a DFT X[ω] de sinais reais é anti-simétrica em torno de πFs .
Para iniciar o estudo das FFTs, consideremos o seguinte sinal periódico (qual o perı́odo
fundamental?)
x(t) = cos(2πt)sen (6πt)
o qual será amostrado com uma taxa Fs . Para conseguir representar adequadamente o sinal
em tempo discreto, essa taxa deve ser superior a duas vezes, no mı́nimo, a frequência máxima
do sinal. Em termos práticos, deve-se considerar uma taxa, ao menos, 10 a 20 vezes maior. O
código a seguir calcula a FFT do sinal x(t) e apresenta os gráficos de magnitude e fase do sinal
transformado.

1
 close all
clear
clc

Ts = 0.01; % intervalo de amostragem

t = 0: Ts :8;
x = cos (2* pi * t ) .* sin (6* pi * t ) ; % sinal no tempo
X = fft ( x ) ; % transformada

N = length ( t ) ; % numero de
pontos
w = 2* pi *(0:( N -1) ) /( N * Ts ) ; % frequencia
angular

figure
subplot (2 ,1 ,1) , plot (w , abs ( X ) ) ; % magnitude
ylabel ( ' | X | ' )
axis tight
subplot (2 ,1 ,2) , plot (w , rad2deg ( angle ( X ) ) ) % fase
ylabel ( ' \ angle X [ deg ] ' )
xlabel ( ' \ omega [ rad / s ] ' )
axis tight
Execute o código e observe o gráfico obtido. Uma primeira observação é de que o gráfico é
anti-simétrico em torno de πFs . Observamos também que X[ω] apresenta dois picos, um em 4π
[rad/s] e outro em 8π [rad/s]. Esses picos confirmam nosso conhecimento de que x[n] consiste
num sinal periódico. De fato, se plotarmos x(t), veremos que é um sinal periódico com perı́odo
0.5 [s].
figure
plot (t , x )
xlabel ( ' t [ s ] ' )
ylabel ( ' x ( t ) ' )
Calculemos em seguida a DFT de um sinal mais complexo, que corresponde à resposta de
um sistema dinâmico linear sujeito a um ruı́do branco de entrada.
G = tf ([10] ,[1 1 10]) ; % funcao de transferencia do sistema
u =0.2* randn (N ,1) ; % ruido de entrada % N definido
anteriormente
y = lsim (G ,u ,t ' ) ; % resposta do sistema % t definido
anteriormente

figure
plot (t ,u ,t , y )
legend ( ' u ( t ) ' , ' y ( t ) ' )
xlabel ( ' t [ s ] ' )
Em seguida é calculado o espectro de y(t) utilizando o comando fft(·) e os gráficos de
magnitude e fase são traçados.
Y = fft ( y ) ;
figure
subplot (2 ,1 ,1)

2
 semilogx (w , 20* log10 ( abs ( Y ) ) ) ;
ylabel ( ' | Y | [ dB ] ' )
axis tight
subplot (2 ,1 ,2)
semilogx (w , unwrap ( rad2deg ( unwrap ( angle ( Y ) ) ) ) ) ;
ylabel ( ' \ angle Y [ deg ] ' )
xlabel ( ' \ omega [ rad / s ] ' )
axis tight
Neste caso o comando unwrap é utilizado para evitar saltos no diagrama de fase. Sem esse
recurso, o diagrama de fase pode ter uma aparência ruidosa, dado que saltos de múltiplos de
2π são possı́veis quando calculamos a fase pelo comando angle.
Alguns pontos significativos dos gráficos traçados são importantes de destacar:

• via de regra, quando as respostas de sistemas dinâmicos são analisadas, utiliza-se escala
logarı́tmica no eixo das abscissas para cobrir uma maior faixa de frequências;

• um outro motivo para usar a escala logarı́tmica é que o comportamento da resposta tende
a variar pouco para valores de frequência próximos (trace os gráficos utilizando o comando
plot e analise o resultado);

• a máxima frequência representável nos gráficos é πFs [rad/s] e está intimamente ligada
com a taxa de amostragem utilizada;

• a amplitude do gráfico de magnitude é dada em decibéis [dB], cuja definição é

YdB [ω] = 20 log(|Y [ω]|).

Tal medida é adotada para ter uma representação mais significativa dos valores de am-
plitude, evidenciando também valores muito próximos de zero.
Além dos pontos destacados, vale frisar que sistemas dinâmicos são, geralmente, passa-
baixas, ou seja, sinais de entrada com frequências baixas aparecem na saı́da do sistema ao
passo que sinais com frequências altas são atenuados na saı́da.
Os gráficos traçados para a resposta do sistema são comumente denominados Diagramas de
Bode e constituem uma importante ferramenta de análise para sistemas lineares. Esse assunto
é abordado com mais detalhes na Seção 3.

2.1 Diferença entre a DFT e a Transformada de Fourier

Para compreender um pouco melhor as diferenças entre a DFT e a FT (Fourier Transform),
tentemos realizar uma convolução usando fft. Definiremos dois vetores e realizaremos a con-
volução por dois métodos distintos: no primeiro usaremos o comando conv; no segundo vamos
multiplicar as DFTs dos dois sinais e depois tomar a transformada inversa.
a =[1 2 3 4 5];
b =[5 1 2 3 4];
c = conv (a , b ) % convolucao no tempo

A = fft ( a ) ;
B = fft ( b ) ;

C = ifft ( A .* B ) % DFT inversa

3
 Executando os comandos listados podemos ver que os resultados são discrepantes. A expli-
cação para isso é que, enquanto o primeiro método consiste em uma convolução linear, o segundo
consiste numa convolução circular, o que acontece devido ao fato de que a DFT assume que o
sinal discreto no tempo é periódico.
Não obstante, é possı́vel obter a convolução linear a partir da DFT. Para isso, é utilizado
um procedimento conhecido como zero-padding, que consiste em preencher os sinais a e b com
zeros de forma que a convolução circular dê o mesmo resultado que a convolução linear.
la = length ( a ) ;
lb = length ( b ) ;

a =[ a zeros (1 , lb -1) ] % zero - padding

b =[ b zeros (1 , la -1) ] % zero - padding

A = fft ( a ) ;
B = fft ( b ) ;

C = ifft ( A .* B )
Execute o arquivo e observe que a convolução linear foi obtida por ambos os métodos. Usar
a FFT desta forma para realizar convoluções é uma prática comum em filtragem de sinais.
Um detalhe digno de nota sobre a implementação da FFT é que os algoritmos são mais
rápidos quando o número de pontos é uma potência de 2. Por isso, é uma prática comum
completar o sinal com zeros até que o número de pontos atinja a próxima potência de 2.

3 Resposta em Frequência e Diagramas de Bode

O termo “resposta em frequência” se refere à resposta em regime permanente de um sistema
a uma entrada senoidal. Os métodos de resposta em frequência foram desenvolvidos entre às
décadas de 1930 e 1940 por Nyquist, Bode, Nichols e outros. Esses métodos são muito impor-
tantes nas teorias de controle clássico e controle robusto, provendo informações fundamentais
tanto para análise do comportamento do sistema em uma faixa larga de valores frequenciais
quanto para o projeto de controladores.
Uma vantagem dos métodos de resposta em frequência é que seus testes são, em geral,
simples e podem ser realizados com exatidão com auxı́lio de um gerador de sinais senoidais.
Variamos a frequência do sinal de entrada e estudamos a resposta resultante. Além disso, esses
métodos podem ser aplicados diretamente a dados experimentais, isto é, dados obtidos a partir
de medições de sistemas fı́sicos (como fizemos anteriormente).
Funções de transferência senoidais são funções complexas da variável ω. Estas funções são
caracterizadas por um módulo e um ângulo de fase; a frequência ω é um parâmetro indepen-
dente. Dentre os métodos existentes para representação gráfica de uma função de transferência
senoidal, focaremos o estudo nos Diagramas de Bode, que correspondem a dois gráficos, gráfico
de módulo e gráfico de fase, os quais devem ser analisados em conjunto.
Como mencionado os diagramas de Bode podem ser obtidos a partir de dados experimentais,
de forma similar ao realizado na Seção 2. Não obstante, é também possı́vel fazer o traçado dos
diagramas de Bode a partir da função de transferência do sistema em estudo. Iniciaremos nosso
estudo traçando os diagramas de Bode referentes a funções de transferência conhecidas e, em
seguida, veremos como obter funções de transferência aproximadas, para sistemas de primeira
e segunda ordens, a partir de dados experimentais.

4
3.1 O comando bode
Para uma função de transferência genérica a tempo contı́nuo

num(s) βm sm + βm−1 sm−1 + · · · + β1 s + β0

G(s) = =
den(s) αn sn + αn−1 sn−1 + · · · + α1 s + α0

com m, n ∈ N, os seus diagramas de Bode podem ser obtidos utilizando o comando bode() de
três formas distintas.
bode ( G )
bode ( num , den )
bode ( num , den , w )

Observação 1. Apesar dos valores de m e n poderem ser quaisquer naturais, usualmente assume-
se m ≤ n para que a função de transferência seja própria. Funções impróprias acarretam
no aparecimento de derivadas sucessivas da função impulso no tempo, não sendo, portanto,
implementáveis.
A primeira forma do comando bode requer, inicialmente, a definição da função de trans-
ferência G(s) por meio do comando tf(); a segunda maneira assume que são passados dois
vetores correspondentes aos coeficientes do numerador (num) e do denominador (den), do termo
de mais alta ordem para o termo de menor ordem; por fim, a última forma, considera que tam-
bém é passado um vetor de frequências (w), nas quais a função de transferência será analisada.
Como exemplo, vamos executar o seguinte código
G = tf ([10] ,[1 1 10]) ; % funcao de transferencia do sistema

figure
bode ( G )
Veja que os diagramas de módulo e fase são apresentados em uma mesma janela. O mesmo
resultado seria obtido utilizando os outros dois comandos, passando um vetor w conveniente no
terceiro caso. Esse vetor pode ser gerado por meio do comando logspace, pois os diagramas
de Bode são traçados em escala logarı́tmica,
w = logspace ( d1 , d2 , N )
Omitindo o terceiro argumento, o comando logspace(d1, d2) gera 50 pontos igualmente es-
paços em escala logarı́tmica entre as frequências 10d1 e 10d2 . O argumento N, por sua vez, serve
para especificar o número desejado de pontos entre as frequências especificadas. Por exemplo,
rode os seguintes comandos e compare os resultados
w = logspace ( -1 , 2)
w = logspace ( -1 , 2 , 150)
Em alguns casos pode ser interessante obter apenas os valores de magnitude e fase para
diferentes valores de frequência e não traçar os diagramas de Bode. Esse comportamento pode
ser obtido do comando bode especificando parâmetros de retorno
[ mag , phase , w ] = bode (.)
Nessa situação, os diagramas não são apresentados e são retornados os valores de magnitude
(mag), em escala linear, e fase (phase), em graus, avaliados nos valores de frequência w. Isso
pode ser conveniente, por exemplo, quando queremos traçar os diagramas de Bode estimados
sobre os gráficos de valores experimentais colhidos (como veremos a seguir).

5
3.2 Diagramas de Bode a partir de dados experimentais
A partir da saı́da medida de um experimento, e de posse da sequência de entrada que originou
tal saı́da, podemos traçar os diagramas de Bode de modo a estimar a função de transferência
do sistema em estudo.
Como exemplo, considere o script apresentado na sequência.
close all
clear
clc

Ts = 0.01; % intervalo de amostragem

t = 0: Ts :8;
N = length ( t ) ;
w = 2* pi *(0:( N -1) ) /( N * Ts ) ; % frequencia angular

G = tf ([10] ,[1 1 10]) ; % funcao de transferencia do sistema

u = chirp (t ,0 , t ( end ) ,50) ' ; % sinal de entrada entra 0 Hz e 50 Hz
y = lsim (G ,u ,t ' ) ; % resposta do sistema

Y = fft ( y ) ;
U = fft ( u ) ;

G2 = Y ./ U ; % funcao de transferencia estimada

[ mag , pha , wg ] = bode ( G ) ;

figure
subplot (2 ,1 ,1)
semilogx ( w (1: end /2) ,20* log10 ( abs ( G2 (1: end /2) ) ) , ' r ' )
hold on
semilogx ( wg ,20* log10 ( squeeze ( mag ) ) , ' b ' )
ylabel ( ' | G | ' )
legend ( ' estimado ' , ' conhecido ' )

subplot (2 ,1 ,2)
semilogx ( w (1: end /2) , rad2deg ( unwrap ( angle ( G2 (1: end /2) ) ) ) , ' r ' )
hold on
semilogx ( wg , squeeze ( pha ) , ' b ' )
ylabel ( ' \ angle G ' )
xlabel ( ' \ omega [ rad / s ] ' )
Execute o arquivo e observe os resultados. Note que os diagramas de Bode obtidos a partir
dos sinais de entrada u e saı́da y estão muito próximos dos diagramas obtidos a partir do
comando bode. Isso ocorre pois estamos trabalhando em um cenário idealizado no qual não há
influência de ruı́dos externos. Para tratar uma situação mais realı́stica, considere que o sinal
de entrada u está corrompido por um ruı́do branco de amplitude 0.01, ou seja,
uc = u + 0.01* randn (N ,1) ;
Tomando uc como o novo sinal de entrada, repita o script acima e analise os resultados. É
interessante executar o arquivo várias vezes, pois, como a entrada uc é sempre diferente, você
notará que algumas vezes a aproximação da função de transferência é melhor do que outras.

6
4 Exercı́cios
1. Gere o sinal x1 [n] = 0.9δ[n], 0 ≤ n ≤ 19. Plote o gráfico de magnitude e fase da FFT
calculada. Os resultados estão coerentes com o que você esperava? Comente.

2. Gere agora o sinal x2 [n] = 0.9δ[n − 5], 0 ≤ n ≤ 19. Plote o gráfico de magnitude e fase
da FFT calculada. Comente.

3. Gere um sinal do seu interesse e plote a magnitude e fase da FFT.

4. Construa e comente o diagrama de Bode de módulo e fase de:
100(s + 5)
a) G(s) =
s2
+ 10s + 100
−100(s + 5)
b) G(s) = 2
s + 10s + 100
25
c) G(s) = 2
s + 4s + 25
9(s2 + 0, 2s + 1)
d) G(s) = , no intervalo 0, 01 ≤ ω ≤ 1000
s(s2 + 1, 2s + 9
5. Neste exercı́cio o intuito é identificar dois sistemas lineares e determinar suas respectivas
funções de transferência a partir de sinais de entrada e saı́da obtidos experimentalmente.
É escolhido um sinal de entrada rico em frequências para excitar o sistema na maior banda
de frequências possı́vel. A saı́da está sujeita a um ruı́do de medida da ordem de −40dB.
O sinal de entrada foi passado por uma função janela para suavizar os efeitos do número
finito de amostras. Os dados para dois sistemas diferentes são fornecidos nos arquivos
[Link] e [Link]. Para cada um dos arquivos, siga as instruções abaixo e apresente
como solução a função de transferência estimada para o sistema.

Instruções:
• Mova os arquivos [Link] e [Link] para seu diretório de trabalho no Matlab.
Cada arquivo contém três variáveis. A variável u representa uma amostra de um
sinal de entrada aplicado a um sistema linear de tempo contı́nuo. A variável y
representa a respectiva saı́da do sistema e a variável time representa o tempo em
que foram tomadas as amostras.
• No Matlab, use o comando load para carregar as variáveis deste arquivo no seu
workspace (ex.: load exp1).
• Calculando as FFTs de u e y, obtenha a função de tranferência H(jω) = Y (jω)/U (jω)
dividindo as duas FFTs.
• Plote o diagrama de Bode da função de transferência obtida com o eixo das frequên-
cias em [rad/s] (ω = 2πf ) devidamente ajustado.
• Sabemos que o sistema pode ser de 1a ou 2a ordem e que não há zeros finitos. A partir
do diagrama de Bode, identifique a ordem do sistema. Lembre-se de que sistemas
de 1a ordem sem zeros finitos apresentam um decaimento de 20 dB por década nas
altas frequências e que sistemas de 2a ordem possuem um decaimento de 40 dB por
década.
• A partir do diagrama, identifique a função de transferência do sistema

7
 – Em se tratando de um sistema de 1a ordem, identifique K e a de forma que
a
H(s) = K
s+a
– Em se tratando de um sistema de 2a ordem, identifique K, ξ e ωn de forma que

ωn2
H(s) = K 2 .
s + 2ξωn s + ωn2

Para tanto, use o fato de que o pico de ressonância Mp e a frequência de resso-

nância ωr obedecem às seguintes relações:

| ln Mp |
ξ=p
π 2 + (ln Mp )2
p
ωr = ωn 1 − 2ξ 2

em que Mp = maxω |H(jω)|/|H(j0)| e ωr = arg maxω |H(jω)|/|H(j0)|.

• Use o comando lsim para obter a resposta à entrada u do sistema com os parâmetros
identificados no item anterior. Plote no mesmo gráfico a saı́da obtida experimental-
mente e a obtida por simulação. Comente.