Laboratório de Sistemas e Sinais
Análise Espectral
Luı́s Caldas de Oliveira
Abril 2009
O objectivo deste trabalho de laboratório é o de ensinar a analisar sinais no domı́nio da frequência. Utilizaremos
dois métodos. O primeiro consiste em representar graficamente os coeficientes da série discreta de Fourier de sinais
de duração finita. O segundo método representa os coeficientes de segmentos de duração limitada de sinais que
variam no tempo, criando o que é chamado de espectrograma.
1 Introdução
Um sinal em tempo discreto de duração finita com p amostras pode ter a seguinte expansão em série de Fourier:
p−1
x(n) = A0 +
2
X
Ak cos(kω0 n + φk )
(1)
Ak cos(kω0 n + φk )
(2)
k=1
se p for um número ı́mpar e
p
x(n) = A0 +
2
X
k=1
se p for par.
Um sinal de duração finita pode ser considerado como um ciclo de um sinal periódico com frequência fundamental ω0 = 2π/p em radianos por amostra, ou f = 1/p em Hertz. Neste laboratório assumiremos que p é sempre
par e faremos a representação gráfica de cada uma das componentes em frequência |A0 |, . . . , |A p/2 | para diversos
sinais com o objectivo de compreender o significado desses coeficientes.
De notar que cada |Ak | dá a amplitude da componente sinusoidal do sinal à frequência kω0 = k2π/p, que
tem como unidades radianos por amostra. Para interpretar estes coeficientes poderá ser conveniente converter esta
unidade para Hertz. Se a frequência de amostragem for f s amostras por segundo, então a conversão poderá ser
feita através de:
(k2π/p)[radianos/amostra] f s [amostras/segundo] k f s
=
[ciclos/segundo]
2π[radianos/ciclo]
p
Assim, cada |Ak | dá a amplitude da componente sinusoidal com frequência k f s /p Hz.
Note que o Matlab não tem nenhuma função pré-definida para calcular os coeficientes da série de Fourier,
tendo no entanto uma função que calcula a transformada rápida de Fourier, chamada fft. Esta função pode ser
usada para calcular os coeficientes da série de Fourier através da seguinte função serieFourier:
function [amplitude, fase] = serieFourier(x)
% SERIEFOURIER - Retorna a amplitude e a fase de cada componente
% sinusoidal da expansão em série de Fourier do sinal dado como
% argumento, que é interpretado como um ciclo de um sinal
% periódico. Assume-se que o argumento tem um número de amostras p que
% é par. O primeiro valor de retorno é um vector contendo as
% amplitudes da componentes sinusoidais na expansão em série de
% Fourier com frequências 0, 1/p, 2/p, ... 1/2. O segundo valor de
% retorno é um vector com as fases das componentes sinusoidais. Ambos
% os vectores têm comprimento de (p/2)+1.
1
p = length(x);
f = fft(x)/p;
amplitude(1) = abs(f(1));
upper = p/2;
amplitude(2:upper) = 2*abs(f(2:upper));
amplitude(upper+1) = abs(f(upper+1));
fase(1) = angle(f(1));
fase(2:upper) = angle(f(2:upper));
fase(upper+1) = angle(f(upper+1));
Se se tiver um vector x com comprimento par, pode-se usar a função para obter os coeficientes da DFS:
[A, phi] = serieFourier(x);
Os vectores A e phi contêm a amplitude e a fase de cada coeficiente.
Para representar graficamente as amplitudes dos coeficientes em função da frequência basta fazer:
p = length(x);
frequencias = [0:fs/p:fs/2];
plot(frequencias, A);
xlabel(’frequencia em Hertz’);
ylabel(’amplitude’);
Em que fs terá o valor da frequência de amostragem em amostras por segundo. A linha
frequencias = [0:fs/p:fs/2];
requer uma análise mais cuidada. Produz um vector com o mesmo comprimento de A, ou seja 1 + p/2, em
que p é o comprimento do vector x. Os elementos do vector frequencias são as frequências em Hertz de cada
componente da série de Fourier.
2 Trabalho para os Alunos
1. Considere o sinal produzido da seguinte forma:
t = [0:1/8000:1-1/8000];
x = sin(2*pi*800*t);
Isto corresponde a 8000 amostras de uma sinusóide de 800 Hz amostrada a 8 kHz. Oiça o vector x. Utilize
a função serieFourier descrita anteriormente para representar graficamente a amplitude dos coeficientes
da série de Fourier de x.
2. O sinal da alı́nea anterior pode ser visto como a amostragem da sinusóide contı́nua:
x(t) = sin(2π800t)
Repare que a frequência angular da sinusóide é a derivada em ordem ao tempo do argumento da função seno:
ω=
d
2π800t = 2π800
dt
Considere agora o sinal
y(t) = sin(2π800t2 )
A este sinal dá-se o nome de chirp. A frequência instantânea pode ser obtida pela derivada do argumento
da função seno:
d
ω(t) = 2π800t2 = 4π800t
dt
Num sinal chirp a frequência varia constantemente com o tempo.
Considere a amostragem a 8 kHz de y(t):
2
t = [0:1/8000:1-1/8000];
y = sin(2*pi*800*(t.*t));
Oiça o sinal e represente graficamente os coeficientes da série de Fourier. Que gama de valores toma a
frequência instantânea?
3. Os coeficientes de Fourier que calculou anteriormente, descrevem a gama de frequências do chirp bem, nas
não a sua dinâmica. Represente graficamente os coeficientes da série de Fourier do sinal z dado por:
z = y(8000:-1:1)
Oiça o sinal. Compare o som de z com o de y e compare os gráficos dos coeficientes de Fourier.
4. O sinal chirp tem uma representação em frequência que varia com o tempo. Mais precisamente, existem
certas propriedades do sinal que mudam suficientemente devagar para o nosso ouvido as entender como uma
variação na composição em frequência do sinal em vez de o considerar como pertencente ao próprio sinal
(como o timbre ou conteúdo tonal). Note que o nosso ouvido não é sensı́vel a frequências abaixo dos 30
Hz. Em vez disso, o nosso cérebro entende essas variações como variações na natureza do som e não como
conteúdo no domı́nio da frequência. Os métodos de análise de Fourier usados anteriormente não reflectem
esse fenómeno psico-acústico.
A série de Fourier localizada procura resolver este problema. O sinal chirp tem 8000 amostras num segundo, mas como não ouvimos variações abaixo dos 30 Hz como conteúdo na frequência, pode fazer sentido
re-analisar o sinal ao ritmo de 30 vezes por segundo. Isto pode ser feito com a seguinte função:
function espectrogramaCascata(s, fs, amostrasespectro, numdeespectros)
% ESPECTROGRAMACASCATA - Faz o gráfico 3-D do espectrograma do sinal
% s.
%
% Argumentos:
% s - o sinal.
% fs - frequência de amostragem em amostras por segundo.
% amostrasespectro - o número de amostras usadas para calcular cada
% espectro.
% numdeespectros - número de espectros a calcular.
frequencias = [0:fs/amostrasespectro:fs/2];
offset = floor((length(s)-amostrasespectro)/numdeespectros);
for i=0:(numdeespectros-1)
start = i*offset;
[A, phi] = serieFourier(s((1+start):(start+amostrasespectro)));
amplitude(:,(i+1)) = A’;
end
waterfall(frequencias, 0:(numdeespectros-1), amplitude’);
xlabel(’frequencia’);
ylabel(’tempo’);
zlabel(’amplitude’);
Esta função pode ser chamada do seguinte modo:
t = [0:1/8000:1-1/8000];
y = sin(2*pi*800*(t.*t));
espectrogramaCascata(y, 8000, 400, 30);
Que produz o gráfico da figura 4. O gráfico mostra 30 conjuntos distintos de coeficientes de Fourier, cada
um calculado com 400 das 8000 amostras disponı́veis. Explique como é que este gráfico descreve o que
ouviu. Crie um gráfico semelhante para o chirp invertido z.
3
0.6
amplitude
0.5
0.4
0.3
0.2
0.1
0
30
4000
20
3000
2000
10
1000
tempo
0
0
frequencia
Figura 1: Representação da série de Fourier localizada de um sinal chirp
5. A figura 4 é fácil de interpretar graças à estrutura relativamente simples do sinal chirp. Sinais mais interessantes são mais difı́ceis de analisar desta forma. Uma forma de visualização alternativa do conteúdo em
frequência é o espectrograma. Um espectrograma é um gráfico como o da figura 4, mas visto de cima. A
altura de cada ponto é representada por uma cor diferente (ou intensidade numa imagem a preto-e-branco).
No Matlab existe uma função pré-definida para gerar um espectrograma:
specgram(y,512,8000);
Isto resulta na imagem apresentada na figura 5. Nessa imagem utilizou-se o mapa de cores por omissão
(jet). Pode experimentar com outros mapas de cores usando o comando colormap. Um particularmente
útil é:
colormap(hot)
Crie uma imagem semelhante para o chirp invertido z. Determine a gama de variação da frequência instantânea.
6. Junto a este relatório encontra alguns ficheiros de áudio. Use os seguintes comandos para os ouvir e visualizar:
[y,fs] = wavread(’audio1.wav’);
soundsc(y,fs)
subplot(2,1,1); specgram(y,1024,fs,[],900)
subplot(2,1,2); plot(y)
4
4000
3500
3000
Frequency
2500
2000
1500
1000
500
0
0
0.1
0.2
0.3
0.4
0.5
Time
0.6
0.7
0.8
0.9
Figura 2: Representação da série de Fourier localizada de um sinal chirp
Interprete os resultados
7. Para o sinal chirp usado anteriormente:
t = [0:1/8000:1-1/8000];
y = sin(2*pi*800*(t.*t));
produza os coeficientes da série de Fourier usando a função serieFourier. Escreva uma função Matlab
que use a equação 2 para reconstruir o sinal original a partir dos coeficientes. A sua função Matlab deverá
começar da seguinte forma:
function x = reconstroi(amplitude, fase)
% RECONSTROI - Dado um vector de amplitudes e um vector de fases,
% constroi um sinal que tem estes valores como coeficientes da série
% de Fourier. Assume-se que os argumentos têm comprimento ı́mpar,
% p/2+1, e que o vector de retorno tem comprimento p.
E deverá realizar a seguinte equação:
∀n ∈ {m ∈ š|m < p}, x(n) =
p/2
X
Ak cos(2π fk n + φk )
k=0
em que Ak e φk são respectivamente a amplitude e a fase dos coeficientes da série de Fourier. Tenha em
atenção que os ı́ndices dos vectores em Matlab começam em 1.
Note que esta função requere um número elevado de operações. Se o seu computador não for suficientemente
potente, construa os coeficientes de Fourier para as primeira 1000 amostras em vez das 8000 e reconstrua
5
o sinal a partir desses coeficientes. Para verificar que a reconstrução funciona, subtraia o sinal reconstruı́do
do sinal original e examine a diferença. A diferença poderá não ser exactamente zero, mas deverá ser muito
pequena quando comparada com o sinal original. Desenhe o gráfico do sinal de diferença.
8. Iremos agora estudar sinais de batimento que correspondem à combinação de sinais sinusoidais com frequências
próximas. Comece por usar as relações de Euler para mostrar que:
2 cos(ωc t) cos(ω∆ t) = cos((ωc + ω∆ )t) + cos((ωc − ω∆ )t)
em que ωc , ω∆ , t ∈ ’.
Esta identidade significa que a multiplicação de dois sinais sinusoidais com frequências ωc e ω∆ é igual à
soma de duas sinusóides com frequências ωc + ω∆ e ωc − ω∆ .
9. Construa um sinal com a soma de duas sinusóides de frequências 790 e 810 Hz, amostradas à frequência
de 8 kHz e com a duração de 1 segundo. Oiça o sinal resultante e descreva o que ouve. Desenhe o gráfico
das primeiras 800 amostras. Mostre como é que o gráfico ilustra o que ouviu e utilize a identidade da alı́nea
anterior para explicar o gráfico.
10. Qual é o perı́odo do sinal da alı́nea anterior? Qual é a frequência fundamental da sua expansão em série de
Fourier? Apresente o gráfico da amplitude dos seus coeficientes de Fourier usando a função serieFourier.
Desenhe o espectrograma usando specgram. Escolha cuidadosamente os parâmetros de specgram para a
imagem ser mais clara. Qual dos dois gráficos representa melhor o que ouviu?
6