Conceitos de Sinais e Sistemas
Mestrado em Ciências da Fala e da Audição
António Teixeira
13 Janeiro 2007
AT 2007
1
Aula
•
•
•
Análise LPC
Análise Cepstral
Obtenção de F0 e
Formante
•
MATLAB
– lpc()
– rceps()
– xcorr()
AT 2007
2
Análise LPC
Uma introdução
AT 2007
3
• A análise de Fourier não é a única forma de
determinar o espectro de um sinal
• Uma técnica muito utilizada na área do
processamento de voz e Fonética envolve determinar
os chamados coeficientes de predição linear (Linear
Predictive Coding Coefficients)
• Este procedimento é conhecido por análise LPC
– é um processo mais complexo que a DFT
– mas é possível compreender os princípios sem entrar nos
detalhes matemáticos mais complexos
AT 2007
4
• Segundo a teoria-fonte filtro produz-se um som pela
passagem de uma excitação por um filtro
“Entrada nula”  cordas vocais  tracto vocal  rad  voz
isto é
“Entrada nula”  sistema (cordas vocais + tracto + rad) voz
• A noção base da análise LPC baseia-se no processo
inverso
– voz  sistema inverso (filtro LPC)  “saída zero”
a saída será zero se o filtro LPC for exactamente o inverso do sistema
• No primeiro caso temos “síntese” ou produção, no
segundo “análise”
AT 2007
5
• Na abordagem LPC as características
espectrais da fonte glotal e radiação são
incluídas conjuntamente com as relativas ao
tracto vocal num mesmo filtro
AT 2007
6
• A análise LPC envolve determinar um filtro
cujas características em termos de resposta em
frequência seja o inverso do espectro do sinal
de voz
• Como já vimos os filtros digitais são definidos
por um conjunto de coeficientes
– lembra-se dos vectores usados no comando filter ?
AT 2007
7
• Também é possível usar um conjunto de
coeficientes para “prever” o valor de uma
amostra do sinal com base em amostras
anteriores
– y[n]= função de y[n-1], y[n-2], y[n-3] ...
• O algoritmo LPC faz uma previsão desta
forma usando um número reduzido de pontos
anteriores, multiplicando cada por um
coeficiente
AT 2007
8
• O princípio básico da análise LPC é a de que
uma amostra pode ser considerado como
simplesmente a soma de um número de
amostras anteriores, cada multiplicada por um
número adequado
– os números são denominados coeficientes de
predição linear
– y[n]= a1 y[n-1] + a2 y[n-2] ....
AT 2007
9
• Para um sinal y, a análise LPC calcula os
coeficientes a[1] ... a[p] tais que
y[n]= a[1] y[n-1] +a[2] y[n-2]+ ... +a[p] y[n-p] +erro
AT 2007
10
Relação com a produção de voz
• A produção pode ser descrito pela equação de convolução
y*a=b*x
onde x é ma fonte, a e b coeficientes.
• Para sons sem anti-ressonâncias (exemplo: as vogais) b=1 e a
equação reduz-se
– y*a=x
y[n] + a[1] y[n-1] +a[2] y[n-2] + ... +a[k]y[n-k]=x[n]
ou rearranjando
y[n] = -a[1] y[n-1]-a[2] y[n-2]- ... -a[k]y[n-k]+x[n]
• Ou seja, num modelo sem anti-ressonâncias, a amostra actual é
igual a uma combinação linear de amostras anteriores da saída
mais a entrada
• O modelo proposto antes é, portanto, adequado
AT 2007
11
Um exemplo de como determinar os
coeficientes
• Consideremos um sinal
• tomemos 12 amostras (uma janela)
– designados por s1, s2, ..., s12
• a estimativa para s5 considerando as 4
amostras anteriores
– s^5= a1 x s4 + a2 x s3 + a3 x s2 + a4 x s1
– num caso concreto teríamos algo como
– s^5 = -42 a1 + 17 a2 + 5 a3 + 90 a4
AT 2007
12
• continuando...
–
–
–
–
s^6= - 71 a1 – 42 a2 + 17 a3 + 50 a4
s^7= ....
s^8= ....
s^9= ....
• Se cada amostra fosse correctamente predicta,
não haveria diferença entre s e s^
0=sn – s^n
AT 2007
13
• isto é
s6 – s^6 = 0
= -40 – (-71 a1 –42 +17 a3 + 50 a4)
s7-s^7=0
= -4 –(-74 a1 –54 a2 + 16 a3 + 97 a4)
s8- s^8=0
= 22 – (-40 a1 – 79 a2 – 54 a3 + 16 a4)
s9-s^9=0
=49 – (-4 a1 –59 a2 – 79 a3 – 54 a4)
temos um sistema de 4 equações com 4 incógnitas
facilmente se obtem a1,a2,a3 e a4
neste caso: a1=0.5, a2=-0.6, a3=0.4 e a4=-0.7
AT 2007
14
generalizando
• No entanto pretendemos obter os coeficientes
que sejam adequados não apenas a este
conjunto restrito e específico de pontos mas
para qualquer amostra
– como vimos resolvendo a equação anterior
obtemos coeficientes adequados para s6 a s9
– mas darão um erro se aplicados a outras amostras
– o erro para cada ponto é designado por en
• en=(s^n – sn)2
– usa-se o quadrado para que seja sempre positivo
AT 2007
15
• teremos:
e6=(s^6 – s6)2
(a1 s5 + a2 s4 + a3 s3 + a4 s2 – s6)2
e7=(s^7 – s7)2
(a1 s6 + a2 s6 + a3 s4 + a4 s3 – s7)2
...
e12=(s^12 – s12)2
(a1 s11 + a2 s10 + a3 s9 + a4 s8 – s7)2
• o algoritmo resolve este conjunto de equações
tentando minimizar o erro
– Usando o Matlab obtêm-se rapidamente usando o
comando lpc()
– Os coeficientes são uma forma eficiente de descrever o
sinal de voz
AT 2007
16
Matlab – lpc()
• A = LPC(X,N) finds the coefficients, A=[ 1 A(2) ...
A(N+1) ], of an Nth order forward linear predictor
– Xp(n) = -A(2)*X(n-1) - A(3)*X(n-2) - ... - A(N+1)*X(n-N)
– such that the sum of the squares of the errors
err(n) = X(n) - Xp(n) is minimized.
• [A,E] = LPC(X,N) returns the variance (power) of the
prediction error.
• LPC uses the Levinson-Durbin recursion to solve the
normal equations that arise from the least-squares
formulation.
– This computation of the linear prediction coefficients is
often referred to as the autocorrelation method.
AT 2007
17
Métodos de obtenção dos coef.
• Existem várias formas de obter os coeficientes
• Sem entrar nos detalhes, refiram-se:
– método da autocorrelação
– método da covariância
– método da “lattice”
– Consultar livros como Rabiner & Schafer 1978
para os detalhes
AT 2007
18
Demo
lpclearn
AT 2007
19
Exemplo
• Consideremos um segmento de sinal (de uma vogal)
• ak=lpc(frame,12)
• resultado
1.0000 -2.3994 2.0545 -0.5626
-0.1950 0.0965 -0.0118 -0.1951
0.7283 -0.5366 -0.4562 0.8068 -0.3069
AT 2007
20
Resposta em frequência
• Como após a determinação dos coeficientes
temos um filtro (em que a saída depende de
valores da saída em instantes anteriores)
podemos obter a sua resposta em frequência
AT 2007
21
• Como o filtro obtido será o inverso do filtro de
produção, as características espectrais do
segmento analisado serão:
AT 2007
22
Análise na frequência com LPC
• Tendo os coeficientes a1...ak facilmente se
obtém o espectro
– Exemplo:
• Material analisado: pequeno segmento de uma vogal
a0= 1.0000 a1= -0.1155
a2= -0.4197 a3= 0.1063
sinal no tempo
0.1
0.05
0
a4= 0.2854
a5= 878
0.6263
Raiz1
Raiz2
a6= -0.2841
a7=1420
-0.2171
Raiz3 2729
a8= -0.0904 a9 = 0.2207
Raiz4 3446
a10= 0.2150 a11= -0.2755
-0.05
-0.1
0
50
100
150
200
250
300
espectro
30
20
10
0
-10
-20
a12= -0.3004
AT 2007
0
500
1000
1500
2000
2500
3000
3500
4000
4500
5000
23
LPC vs FFT para obter espectro
Espectro usando FFT e LPC
30
20
10
0
-10
-20
-30
-40
-50
0
AT 2007
500
1000
1500
2000
2500
3000
3500
4000
4500
5000
24
• A análise LPC separa os componentes
relativos à fonte e ao filtro
• É importante para a determinação da
frequência fundamental e as formantes
AT 2007
25
Questões práticas
• Ordem a utilizar
– Regra prática
• Frequência de amostragem em kHz + 2
• Exemplo: 10000 Hz => 10+2=12
• Aplicar janelas
• Usar pré-ênfase
– O espectro da fala decai com o aumento da
frequência
– Para facilitar a análise LPC tenta-se corrigir esse
efeito
• y(n)=x(n) – a x(n-1), a~0.98
AT 2007
26
Pré-ênfase
• As primeiras formantes têm maior energia e são
preferencialmente modeladas
– A maior energia deve-se ao efeito combinado da excitação
glotal e da radiação
• Geralmente utiliza-se um filtro de pré-ênfase
– s’(n)=s(n) – a1 s(s-1)
– Tipicamente 0.96<= a1 <= 0.99
• Para reconstruir o sinal deve usar-se o filtro inverso
– s(n)= s’(n) + a1 s(n-1)
AT 2007
27
Leitura adicional
• Capítulo 11 do livro “Elements of Acoustic
Phonetics” de Peter Ladefoged, 2ª ed., University of
Chicago Press.
– SDUA 801.4 17 2ed
• Capítulo 8 do livro “Techniques in Speech Acoustics”
de J. Harrington e S. Cassidy, Kluwer Academic
Press, 1999
• SDUA 800H 664
• Apresenta:
– informação sobre a forma como são calculados os coeficientes LPC
(secção 8.2)
– Obtenção do espectro com base nos coeficientes (sec. 8.4)
AT 2007
28
Exercícios Matlab
• Obter um pequeno segmento (frame) de uma vogal
• Obter os coeficientes com a ajuda do Matlab
• Obter a resposta em frequência do filtro e do filtro
inverso
– comparar com o espectro obtido pela DFT/FTT
• Obter o sinal de erro
• Verificar o efeito de alterar o número de coeficientes
• Repetir o processo para outro tipo de som (fricativa
por exemplo)
AT 2007
29
Análise Cepstral
AT 2007
30
Motivação
• Como o sinal de voz pode ser obtido pela convolução
da excitação glotal com a resposta impulsional do
filtro constituído pelo tracto torna-se necessário
muitas vezes efectuar a operação inversa
(desconvolução)
• A análise cepstral é uma das técnicas que permite
estimar uma separação da fonte do filtro
• Uma das motivações é que os harmónicos da
frequência fundamental podem dificultar a análise das
formantes
– uma muito melhor estimativa das formantes poderia ser
obtida se os harmónicos forem removidos de alguma forma
AT 2007
31
Propriedades importantes
• Para compreender a análise cepstral interessa
perceber como são representados no espectro a fonte
e o filtro
• Uma das propriedades da DFT é que se dois sinais x
(a fonte) e h (o filtro) são convoluidos a sua DFT é
igual ao produto da DFT de x pela DFT de h
• Quando se representa o espectro em dBs temos uma
escala logaritmica
– logaritmo(a x b) = logaritmo(a) + logaritmo (b)
• Portanto, o espectro em dBs representa a SOMA do
espectro da fonte com o do filtro
– o que nos fornece um caminho para os separar ...
AT 2007
32
• Uma pequena revisão:
– se tivermos um sinal composto por duas
sinusóides, uma variando lentamente no tempo,
outra rapidamente,
• isto é uma de baixa frequência e outra de frequência
elevada
– onde apareceriam as riscas no espectro ?
– A correspondente à baixa frequência apareceria na
parte “baixa” do espectro; a outra na parte “alta”
AT 2007
33
O cepstro
• É esta lógica que está na base da separação das
variações rápidas do espectro devido à fonte das
variações lentas do filtro
• Se considerarmos o espectro como um sinal (no
tempo) e aplicarmos a DFT então a parte devida à
fonte deverá aparecer nas frequências elevadas e a
relativa ao filtro nas frequências baixas
– Na prática não se aplica a DFT para a inversa (IDFT) para
converter da frequência para o tempo
• o caminho inverso da DFT
– Curiosamente, apesar de inversas a DFT e a IDFT resultam
no mesmo efeito de separação
• com a DFT ou IDFT a parte de variação rápida é separada da parte
de variação lenta
AT 2007
34
Cepstro real
Definição:
AT 2007
35
Origem do nome
• Como se trata de um espectro de outro
espectro, os seus “inventores” criaram uma
variação da palavra “spectrum” chegado a
“cepstrum”
– adaptado ao Portugês: “espectro”  “cepstro”
• Já agora:
– filtering  liftering
– frequency  quefrency
AT 2007
36
Em Matlab
0.3
0.2
0.1
x
0
-0.1
-0.2
-0.3
-0.4
0
10
20
30
40
50
60
70
z=rceps(x);
0.5
0
-0.5
-1
z
x
-1.5
-2
AT 2007
0
50
100
150
200
250
37
Espectro suave
• Depois de separados podemos eliminar cada
uma das partes por um processo de filtragem
– Se eliminarmos a parte de frequências mais
elevadas e voltarmos a efectuar uma DFT teremos
um espectro “suave” com “apenas” as
características devidas ao filtro (tracto)
AT 2007
38
Cepstro de sinal periódico
• Se o sinal original é periódico então a fonte
manifesta-se como “picos” espaços da
duração do período fundamental
AT 2007
39
Obtenção da Frequência Fundamental
AT 2007
40
Determinação de F0
• F0 é uma propriedade fundamental dos sons
vozeados
• Estimar F0 é muito mais difícil do que se
possa imaginar !!
– A excitação é apenas quase-periódica
• Alguns Métodos
– Método da autocorrelação
– Método usando predição linear
– Método cepstral
AT 2007
41
Pitch e frequência
• Pitch é a qualidade subjectiva relacionada com
a frequência
– No entanto, outros factores afectam a percepção de
pitch
• por exemplo: o pitch depende em certa medida da
intensidade com que um tom é apresentado ao ouvinte
AT 2007
42
Determinação do pitch pelo método da autocorrelação
1
0.5
0
-0.5
0
2
4
6
8
10
12
14
16
A janela deve conter pelo menos dois
períodos de pitch
AT 2007
18
20
close all;clear all
[x,fs]=wavread('seg4');
t=(1:length(x))/fs*1000;
plot(t,x)
%Defina janela de observaçao de 20ms
N=floor(0.02*fs);
t1=(1:N)/fs*1000;
rx=xcorr(x,N,'coeff');
figure(2)
plot(t1,rx(N+1:2*N))
%determine o maximo da autocorrelaçao para
%desvios superiores a 2ms(500Hz)
N1=floor(0.002*fs);
[x0,imax]=max(rx(N+N1:2*N+1));
imax=imax+N1;
t0=imax/fs*1000;
f0=1/t0*1000;
fprintf(1,'O pitch e´: %6.2f ms\n',t0)
fprintf(1,'A frequencia fundamental e´: %6.1f Hz\n',f0)
O pitch e´: 8.00 ms
A frequencia fundamental e´: 125.0 Hz
43
Determinação da autocorrelação
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
0
50
100
150
• Estimativa biased
200
250
300
350
400
450
500
• Estimativa unbiased
-3
6
12
5
x 10
10
4
8
3
6
2
4
1
2
0
0
-2
-1
-4
-2
-6
-3
AT 2007
0
100
200
300
400
500
600
700
800
900
1000
-8
0
100
200
300
400
500
600
700
800
900
1000
44
• A utilização directa da autocorrelação pode
resultar em múltiplos máximos
– Tornando difícil a decisão
• Um método para tentar resolver o problema é
utilizar “center-clipping”
– Colocando a zero as amostras que se situem abaixo
de um certa percentagem da amplitude máxima
(por exemplo Sondhi usou 30 %)
AT 2007
45
exemplo
autocorrelação
0.08
0.06
0.04
1
0.02
0
-0.02
-0.04
0.5
-0.06
-0.08
0
200
400
600
800
1000
1200
1400
0
Center-clipped 50 %
-0.5
0
2
4
6
8
10
12
14
16
18
20
O pitch e´: 4.60 ms
A frequencia fundamental e´: 217.4 Hz
AT 2007
46
Determinação do pitch pelo cepstrum
AT 2007
47
Determinação do pitch pelo cepstrum
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
0
10
20
30
40
50
60
70
0.5
0
-0.5
-1
-1.5
-2
0
50
100
150
200
250
close all;clear all
[x,fs]=wavread('seg4');
N=length(x);
t=(1:length(x))/fs*1000;
plot(t,x)
z=rceps(x);
figure(2)
plot(z(1:length(x)/2))
N1=0.02*N
[z0,imax]=max(z(N1:N/2));
imax=imax+N1
t0=imax/fs*1000;
f0=1/t0*1000;
fprintf(1,'O pitch e´: %6.2f ms\n',t0)
fprintf(1,'A frequencia fundamental e´: %6.1f
Hz\n',f0)
O pitch e´: 7.88 ms
A frequencia fundamental e´: 127.0 Hz
AT 2007
48
Outro exemplo
• Mesmo sinal usado em center-clipped
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
-1.2
-1.4
-1.6
0
100
200
300
400
500
600
700
O pitch e´: 4.59 ms
A frequencia fundamental e´: 218.1 Hz
AT 2007
49
• O cepstrum contém harmónicos da frequência
fundamental
• Os valores baixos de quefrency representam a
forma do tracto
• Os valores elevados de quefrency representam
a excitação
– E no caso de sinais vozeados a frequência
fundamental
AT 2007
50
AMDF
• AMDF – Average Magnitude Difference
Function
1
N
N
 | s ( n)  s ( n  i ) |
i 0
• Mais rápido, em especial quando se utiliza
aritmética inteira
– Não necessita de multiplicações
AT 2007
51
Determinação do pitch por filtragem inversa
Filtro
passa
baixo
Filtro
inverso
Janela
Autocor
relação
Análise
LPC
-3
0.25
4
0.2
x 10
1
3
0.15
0.8
2
0.6
0.1
1
0.05
0.4
0
0
-1
-0.05
-2
-0.1
0
-0.15
-3
-0.2
-4
-0.25
0.2
0
5
10
15
20
25
30
-5
-0.2
0
2
4
6
8
10
12
14
16
18
20-0.4 0
2
4
6
8
10
12
14
16
18
20
O pitch e´: 8.13 ms
A frequencia fundamental e´: 123.1 Hz
AT 2007
52
Determinação do pitch por filtragem inversa
close all;clear all
[x,fs]=wavread('seg1');
t=(1:length(x))/fs*1000;
plot(t,x)
%filtragem passa baixo
[b,a]=butter(3,0.25);
x=filter(b,a,x);
%Defina janela de observaçao de 20ms
N=floor(0.02*fs);
y=x(1:N).*hamming(N);
t=(1:N)/fs*1000;
%Determine o modelo LPC de
ordem 16
p=16;
a=real(lpc(y,p));
%determinação do residuo por
filtragem
%inversa
e=filter(a,1,y);
AT 2007
figure(2)
plot(t,e)
ry=xcorr(y,N,'coeff');
figure(3)
plot(t,ry(N+1:2*N))
%determine o maximo da
autocorrelaçao para
%desvios superiores a 2ms(500Hz)
N1=floor(0.002*fs);
[x0,imax]=max(ry(N+N1:2*N+1));
imax=imax+N1;
t0=imax/fs*1000;
f0=1/t0*1000;
fprintf(1,'O pitch e´: %6.2f ms\n',t0)
fprintf(1,'A frequencia fundamental e´:
%6.1f Hz\n',f0)
53
Pós-processamento
• Os métodos expostos podem cometer erros
– Produzindo variações bruscas do valor do pitch
que são incorrectas
• Muitas vezes recorre-se a pós-processamento
– Filtro de mediana
•
•
•
•
Filtro de comprimento L (3 ou 5)
entrada L valores de pitch
saída a mediana (L-1)/2 valores abaixo, (L-1)/2 valores acima
Pode usar-se um filtro passa baixo depois do filtro de mediana
– Programação dinâmica
• Algoritmo de optimização
AT 2007
54
F0 usando SFS
AT 2007
55
Formantes
AT 2007
56
Porquê calcular as formantes ?
• As formantes são definidas perceptualmente
• A propriedade física correspondente é a
frequência de ressonância do tracto vocal
• Análise de formantes é útil para posicionar os
fonemas em termos das primeiras 2 ou 3
formantes
– As duas primeiras formantes
identificam/caracterizam bastante bem as vogais
AT 2007
57
Obter valores candidatos
• Procura de picos no espectro
– Designado em Inglês de “peak picking”
• Procura de picos no espectro obtido de análise LPC
– Várias alternativas:
•
•
•
•
Reter os N maiores picos,
Os N picos com menores frequências
Todos os picos
Pontos onde a segunda derivada é mais negativa
• Factorização das raízes do polinómio resultante da
análise LPC
AT 2007
58
Processos habituais
Pré - ênfase
1-0.95 z-1
Janela
Hamming
Cálculo dos ak
Peak picking
Cálculo de |A(ejw )|2
usando FFT
Cálculo das raizes de A(z)
Procura de máximos
Decisão
1
p
1   ak z  k
k 1
AT 2007

Melhorado pela
utilização de
Interpolação
Parabólica
(Boite et al.
P 92)
p
1
p
 (1  z z
k
1
Ck
1
k 1 (1  z k z )

)
Fk, Bk
k 1
59
Cálculo de Fk e Bk
• Uma raiz
zk  k .e
j k
• Próxima do circulo unitário corresponde a uma
formante, com:
k
Fk 
f amostragem
2
Bk  ( f ,3dB )  1 /  .(1   k ) f amostragem
AT 2007
60
Exemplo “seg6”
F1 = 326.40 Hz
close all;clear all
0.15
[x,fs]=wavread('seg6');
0.1
F2 = 1133.75 Hz
0.05
t=(1:length(x))/fs*1000;
0
%Defina janela de observaçao de 20ms
-0.05
F3 = 2824.89 Hz
N=floor(0.02*fs);
-0.1
y=x(1:N).*hamming(N);
-0.15
-0.2
0
10
20
30
40
F4 = 4039.82 Hz
50
60
70
t=(1:N)/fs*1000;
%Determine o modelo LPC de ordem 12
p=12; a=real(lpc(y,p));
% raizes
1
zplane(1,a); rs=roots(a);
Miuk=abs(rs);tetak=angle(rs)
% eliminar metade
Imaginary part
0.5
ind=find(tetak<=0);miuk(ind)=[];tetak(ind)=[];
% Fk
0
fk=tetak/(2*pi)*fs;
[fk,ind]=sort(fk);
-0.5
% mostar resultados
fprintf(1,'F1 = %6.2f Hz\n',fk(1))
-1
fprintf(1,'F2 = %6.2f Hz\n',fk(2))
-1
-0.5
0
Real part
0.5
1
fprintf(1,'F3 = %6.2f Hz\n',fk(3))
fprintf(1,'F4 = %6.2f Hz\n',fk(4))
AT 2007
61
No SFS
F1= 355
F2=1168
F3=2809
AT 2007
62