Cap13-slides

Prévia do material em texto

Análise Espectral 
 Uma aplicação importante dos métodos de processamento digital de 
sinais é na determinação do conteúdo em frequência de um sinal 
contínuo 
 
 Análise espectral: determinação do espectro de energia ou espectro 
de potência do sinal (PSD – Power Spectral Density) 
 
 Dado um sinal analógico 𝑥𝑎(𝑡), limitado em frequência, este é 
convertido para um sinal digital 𝑥 𝑛 , assumindo-se que o efeito de 
aliasing pode ser ignorado e que a precisão do conversor é 
suficiente para que efeitos de quantização sejam desprezíveis 
 
 A DFT (calculada de forma rápida pela FFT) pode ser usada para a 
análise espectral de sinais de comprimento finito 
Análise Espectral 
 Neste caso, obtemos as amplitudes e fases de cada componente 
senoidal do sinal 𝑥 𝑛 , ou seja: 
 
𝑥 𝑛 =
1
𝑁
 𝑋 𝑘 𝑒𝑗
2𝜋
𝑁 𝑘𝑛
𝑁−1
𝑘=0
 
 Se 𝑥 𝑛 for real: 
𝑥 𝑛 =
1
𝑁
𝑋 0 +
2
𝑁
 𝑋 𝑘 cos 
2𝜋
𝑁
𝑘𝑛 + ∠𝑋 𝑘
𝑁/2
𝑘=1
 
 
 Para aumentar a resolução, pode-se fazer uma interpolação 
espectral, acrescentando-se zeros no final da sequência e 
calculando-se uma DFT de tamanho maior do que o seu 
comprimento. 
 
Análise Espectral 
 Exemplo de interpolação: 
 
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
-40
-30
-20
-10
0
10
20
30
40
Frequência digital normalizada
A
m
p
lit
u
d
e
 (
d
B
)
 
 
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
-40
-20
0
20
40
Frequência digital normalizada
A
m
p
lit
u
d
e
 (
d
B
)
 
 
N=32
N=128
n=[0:31]; 
x=sin(2*pi*7*n/32)+0.2*sin(2*pi*13*n/32)+
0.01*randn(1,32); 
X=fft(x); 
Xdb=20*log10(abs(X)); 
subplot(2,1,1); 
plot(n/16,20*log10(abs(X))) 
xlabel('Frequência digital normalizada') 
ylabel('Amplitude (dB)'); 
legend('N=32') 
axis([0 2 -40 40]) 
 
Xi=fft(x,128); 
k=[0:127]; 
Xidb=20*log10(abs(Xi)); 
subplot(2,1,2); 
plot(k/64,20*log10(abs(Xi)),'r') 
xlabel('Frequência digital normalizada') 
ylabel('Amplitude (dB)'); 
legend('N=128') 
axis([0 2 -40 40]) 
 
Análise Espectral 
 Efeito de “vazamento” (leakage): Quando uma componente senoidal 
tem frequência não múltipla de 2𝜋/𝑁, sendo 𝑁 o comprimento da 
sequência, devido a não haver um número inteiro de períodos desta 
componente 
 
 Exemplo: 
 n=[0:31]; 
x=sin(2*pi*4*n/32)+sin(2*pi*6*n/32)+0.05*sin(2*pi*13*n/32); 
X=fft(x); 
subplot(2,1,1); stem(n/16,abs(X)) 
xlabel('Frequência digital normalizada') 
ylabel('Amplitude (dB)'); 
legend('X') 
axis([0 2 -5 20]) 
 
x2=sin(2*pi*4.5*n/32)+sin(2*pi*6.5*n/32)+0.02*sin(2*pi*13*n/32); 
X2=fft(x2); 
subplot(2,1,2); stem(n/16,abs(X2),'r') 
xlabel('Frequência digital normalizada') 
ylabel('Amplitude (dB)'); 
legend('X2') 
axis([0 2 -5 20]) 
Análise Espectral 
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
-5
0
5
10
15
20
Frequência digital normalizada
A
m
p
lit
u
d
e
 (
d
B
)
 
 
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
-5
0
5
10
15
20
Frequência digital normalizada
A
m
p
lit
u
d
e
 (
d
B
)
 
 
X
X2
Análise Espectral 
 Para reduzir o efeito de “vazamento”, multiplica-se o sinal por uma 
janela, que decai suavemente no início e no final do intervalo de 
observação do sinal 
 
 Janelas mais usadas: Hamming, Hanning, Blackmann e Kaiser 
 
 Com o emprego de janelas: 
 
 Haverá perda de resolução para diferenciar componentes de 
frequências próximas, devido ao aumento da largura do lóbulo 
principal da janela (em relação à janela retangular); para 
aumentar resolução é necessário aumentar o comprimento da 
janela (intervalo de observação) 
 
 Aumentará a discriminação de componentes mais fracas na 
presença de componentes mais fortes, devido à diminuição dos 
lóbulos secundários (em relação à janela retangular) 
 
 
 
Análise Espectral 
 Exemplo de vazamento e do uso de janelas: 
 x=sin(2*pi*4.5*n/32)+sin(2*pi*6.5*n/32)+0.02*sin(2*pi*13*n/32); 
 xh=x.*hanning(32); 
 xb=x.*blackman(32); 
 
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
-80
-70
-60
-50
-40
-30
-20
-10
0
10
Frequência digital normalizada
A
m
p
lit
u
d
e
 (
d
B
)
 
 
Retangular
Hanning
Blackman
Análise Espectral de Sinais Não Estacionários 
 Os parâmetros das componentes senoidais (amplitude, fase e/ou 
frequência) do sinal são variantes no tempo 
 
 Sinais de voz, de radar e de sonar são exemplos de sinais não 
estacionários 
 
 Exemplo: sinal chirp 
𝑥 𝑛 = 𝐴𝑐𝑜𝑠 𝜔0𝑛
2 
 
 obtido do sinal contínuo 𝑥(𝑡) = 𝐴𝑐𝑜𝑠 𝜔0𝑡
2 = 𝐴𝑐𝑜𝑠 𝜙(𝑡) . 
 A frequência instantânea de 𝑥(𝑡) é 
𝑑𝜙(𝑡)
𝑑𝑡
= 2𝜔0𝑡 (linear com o 
 tempo). 
 
 
Análise Espectral de Sinais Não Estacionários 
 Podemos analisar trechos de um sinal não estacionário, multiplicando-o 
por uma janela 𝑤(𝑛) e deslocando-o em relação à janela 
 
 À cada subsequência resultante, é aplicada uma DFT, gerando uma 
transformação denominada STFT 
 
 A STFT (Short-Time Fourier Transform) é definida como 
 
𝑋𝑆𝑇𝐹𝑇 𝑘,𝑚 = 𝑥 𝑚𝐿 + 𝑛 𝑤[𝑛]𝑒
−𝑗
2𝜋𝑘𝑛
𝑁
𝑁−1
𝑛=0
, 0 ≤ 𝑘 ≤ 𝑁 − 1 
 
 Função de duas variáveis: 
 𝑘: índice correspondente à frequência 𝜔𝑘 =
2𝜋𝑘
𝑁
 
 𝑚: índice correspondente ao deslocamento da janela de 
observação de 𝑚𝐿 amostras 
Análise Espectral de Sinais Não Estacionários 
 A STFT é visualizada por um gráfico 3D ou com a intensidade em 
cores (espectrograma) 
 
%EXAMPLE: Display the PSD of each segment of a linear chirp. 
 t=0:0.001:2; % 2 secs @ 1kHz sample rate 
 y=chirp(t,100,1,200); % Start @ 100Hz, cross 200Hz at t=1sec 
 spectrogram(y,128,120,128,1E3); % Display the spectrogram 
 title('Linear Chirp: start at 100Hz and cross 200Hz at t=1sec'); 
 
 S=spectrogram(y,128,120,128,1E3); % Display the spectrogram 
 figure,mesh(abs(S)) 
 
0 50 100 150 200 250 300 350 400 450 500
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Frequency (Hz)
Linear Chirp: start at 100Hz and cross 200Hz at t=1sec
T
im
e
0
50
100
150
200
250
0
10
20
30
40
50
60
70
0
5
10
15
20
25
30
35
Deslocamento da janela (em amostras)
Frequência ( 1000/128)
M
a
g
n
it
u
d
e
Análise Espectral de Sinais Não Estacionários 
 Espectrograma de um sinal de voz 
 
 [x,fs]=wavread('female_src_1'); 
 spectrogram(x,512,256,1024,fs); 
 
Análise Espectral de Sinais Aleatórios 
 A análise espectral de sinais aleatórios pode ser feita utilizando-se o 
método do Periodograma ou método de Welch, que estima a 
densidade espectral de potência (PSD) 
 
 Neste método, calculam-se STFTs de trechos do sinal (geralmente 
com sobreposição) e tira-se a média dos seus quadrados, obtendo-se 
o periodograma do sinal 
 
𝑃𝑥 𝑘 =
1
𝑀
 𝑋𝑆𝑇𝐹𝑇 𝑘,𝑚
𝑀−1
𝑚=0
2
, 0 ≤ 𝑘 ≤ 𝑁 − 1 
 
 Em geral, usam-se janelas para melhorar a discriminação de 
componentes de baixa potência 
 
 Funções no Matlab: periodogram, pwelch 
 
 
Análise Espectral de Sinais Ruidosos 
 Exemplo: 𝑥 𝑛 = 5𝑠𝑒𝑛 0,4𝜋𝑛 + 𝑣 𝑛 
 𝑣 𝑛 : ruído branco de variância unitária 
(a) Periodogramas de 50 
realizações com N=64 e 
(b) Periodograma médio 
(c) Periodogramas de 50 
realizações com N=256 e 
(b) Periodograma médio 
Análise Espectral de Sinais Ruidosos 
 Exemplo: 
 
 
 
Fs = 1000; t = 0:1/Fs:.3; 
x = cos(2*pi*t*200)+randn(size(t)); % A cosine of 200Hz plus noise 
periodogram(x,[],'twosided',512,Fs); % The default window is used 
0 100 200 300 400 500 600 700 800 900
-55
-50
-45
-40
-35
-30
-25
-20
-15
-10
Frequency (Hz)
P
o
w
e
r/
fr
e
q
u
e
n
c
y
 (
d
B
/H
z
)
Periodogram Power Spectral Density Estimate
Análise Espectral de Sinais Ruidosos 
 Exemplo: 
 
 
 
Fs = 1000; t = 0:1/Fs:.296; 
 x = cos(2*pi*t*200)+randn(size(t)); % A cosineof 200Hz plus noise 
 pwelch(x,[],[],[],Fs,'twosided'); % Uses default window, overlap & NFFT 
0 100 200 300 400 500 600 700 800 900
-38
-36
-34
-32
-30
-28
-26
-24
-22
-20
-18
Frequency (Hz)
P
o
w
e
r/
fr
e
q
u
e
n
c
y
 (
d
B
/H
z
)
Welch Power Spectral Density Estimate
Transformada Cosseno Discreta (DCT) 
 A DFT de uma sequência real 𝑥[𝑛] é uma sequência complexa 𝑋[𝑘] 
com simetria em torno da amostra 𝑘 =
𝑁
2
 
 
 A DCT é uma transformação ortogonal que representa uma sequência 
real no domínio do tempo, 𝑥 𝑛 , por uma sequência também real no 
domínio da transformada, 𝑋[𝑘] 
 
 Outras transformadas reais são a DST (discrete sine transform) e a 
Haar 
 
 A DCT pode ser implementada de forma eficiente usando FFTs 
 
 É utilizada em compressão de sinais, principalmente de imagens 
(JPEG), pela sua propriedade de compactação de energia 
 
 
Transformada Cosseno Discreta (DCT) 
 Há diferentes tipos de DCT: 
 DCT-I: 
𝑋 𝑘 =
1
2
𝑥 0 + (−1)𝑘𝑥 𝑁 − 1 + 𝑥 𝑛 𝑐𝑜𝑠
𝜋
𝑁 − 1
𝑘𝑛
𝑁−2
𝑛=1
 
 DCT-II: 
𝑋 𝑘 = 𝑥 𝑛 𝑐𝑜𝑠
𝜋
𝑁
𝑛 +
1
2
𝑘
𝑁−1
𝑛=0
 
 DCT-III: 
𝑋 𝑘 =
1
2
𝑥 0 + 𝑥 𝑛 𝑐𝑜𝑠
𝜋
𝑁
𝑛 𝑘 +
1
2
𝑁−1
𝑛=1
 
 DCT-IV: 
𝑋 𝑘 = 𝑥 𝑛 𝑐𝑜𝑠
𝜋
𝑁
𝑛 +
1
2
𝑘 +
1
2
𝑁−1
𝑛=0