[bibshow file=ref.bib]

Vamos imaginar que ao amostrar um sinal qualquer, tenhamos obtido a seguinte sequência:

$[2.0,\ 3.5,\ 3.9,\ 2.8,\ 0.5,\ -2.9],$

e, num exercício de poderes psíquicos (ou matemáticos, preferívelmente), desejássemos estimar sua próxima amostra. Como poderíamos proceder?

– Chamamos um vidente?

– Usamos o Olho de Agamotto?

– Sentamos e choramos?

Bem, os conceitos de processamento digital de sinais e uma boa dose de álgebra linear podem nos dar alternativas mais viáveis para resolvermos esse problema. Antes, para facilitar e padronizar a notação, vamos chamar essa sequência de x[n], de forma que $x[1]=2.0$ , $x[2]=3.5$ , etc.¹ Uma ideia interessante seria a de utilizar uma combinação linear com amostras anteriores. Dessa maneira, poderíamos estimar a próxima amostra desse sinal usando, por exemplo, 3 amostras anteriores da seguinte maneira:

$\hat x[7] = a_1\times(-2.9) + a_2\times(0.5) +a_3\times(2.8).$

Essas constantes ( $a_1$ , $a_2$ e $a_3$ ) devem ser capazes de representar bem qualquer amostra do sinal com base em suas 3 amostras anteriores. Dessa forma, podemos utilizar as amostras que já conhecemos do sinal para determinarmos os valores dessas constantes. Assim, poderíamos estimar as amostras de $x$ da seguinte maneira²:

$\begin{matrix} \hat x[4] &=& a_1 x[3] + a_2 x[2] + a_1 x[1]\\ \hat x[5] &=& a_1 x[4] + a_2 x[3] + a_1 x[2]\\ \hat x[6] &=& a_1 x[5] + a_2 x[4] + a_1 x[3]\\ \vdots & & \vdots \end{matrix}$

de forma que $a_1$ , $a_2$ e $a_3$ sejam estimados para que $\hat x$ seja o máximo possível similar à $x$ .

Agora, generalizando o nosso problema, vamos assumir que temos um sinal $x[n]$ qualquer, com $N$ amostras, e vamos utilizar um número $p$ qualquer de amostras anteriores para estimar a próxima amostra desse sinal. Dessa forma, teremos que $\hat x[n]$ será dado como segue:

(1) $\begin{equation*} \hat x[n]=\sum_{k=1}^p a_k x[n-k]. \end{equation*}$

Assim, podemos definir o erro proporcionado pela estimação como a diferença entre as amostras de $x$ e as amostras de $\hat x$ , i.e.:

(2) $\begin{equation*} e[n]=x[n] - \hat x[n]. \end{equation*}$

Reescrevendo (2) com base em (1), temos que:

(3) $\begin{equation*} x[n]=\sum_{k=1}^p a_k x[n-k] \ + \ e[n]. \end{equation*}$

Para que nosso preditor funcione bem, precisaremos encontrar esses coeficientes $a_k$ , de forma que tenhamos um menor erro na estimação de $x$ . Mas afinal, como podemos fazer isso?

Visualizar esse problema na forma matricial pode nos ajudar! A equação a diferenças em (3) é equivalente à seguinte expressão matricial:

(4) $\begin{equation*} \begin{bmatrix}x[1]\\x[2]\\x[3]\\ \vdots \\ x[N]\end{bmatrix} = \begin{bmatrix}x[0] & x[-1] & \dots & x[-p+1] \\ x[1] & x[0] & \dots & x[-p+2] \\ x[2] & x[1] & \dots & x[-p+3] \\ \vdots & \vdots & \ddots & \vdots \\ x[N-1] & x[N-2] & \dots & x[N-p]\end{bmatrix}\begin{bmatrix}a_1\\a_2\\ \vdots \\ a_p\end{bmatrix} + \begin{bmatrix}e[1]\\e[2]\\e[3]\\ \vdots \\ e[N] \end{bmatrix}. \end{equation*}$

Se definirmos

$X=\begin{bmatrix}x[1]\\x[2]\\x[3]\\ \vdots \\ x[N]\end{bmatrix}, M=\begin{bmatrix}x[0] & x[-1] & \dots & x[-p+1] \\ x[1] & x[0] & \dots & x[-p+2] \\ x[2] & x[1] & \dots & x[-p+3] \\ \vdots & \vdots & \ddots & \vdots \\ x[N-1] & x[N-2] & \dots & x[N-p]\end{bmatrix}, A=\begin{bmatrix}a_1\\a_2\\ \vdots \\ a_p\end{bmatrix}\ \text{e}\ \ E=\begin{bmatrix}e[1]\\e[2]\\e[3]\\ \vdots \\ e[N] \end{bmatrix},$

podemos reescrever (4) da seguinte maneira:

(5) $\begin{equation*} X=M A + E. \end{equation*}$

Para uma estimação adequada, os valores de $A$ devem minimizar o erro quadrático $\left(\sum {e[k]}^2\right)$ . Para que isso ocorra, de acordo com o método dos mínimos quadrados, temos que $A$ deverá ser estimado da seguinte maneira³:

(6) $\begin{equation*} A = M^\dagger X, \end{equation*}$

(7) $\begin{equation*} A = (M^TM)^{-1}M^T X. \end{equation*}$

Em caso de dúvidas do porquê dessa estimativa para A, clique aqui para mais detalhes

Partindo de (5) temos que:

$E = X - MA.$

Assim, temos que o erro quadrático (Vamos chamá-lo de $S$ ) pode ser calculado como segue (obs: $E^TE=\sum {e[k]^2}$ ):

$\begin{align*} S=E^TE&={(X-MA)}^T(X-MA)\\ &=X^TX-X^TMA-A^TM^TX+A^TM^TMA \end{align*}$

$S=X^TX-2A^TM^TX+A^TM^TMA.$

Como queremos minimizar $S$ , queremos encontrar algum $A$ que satisfaça $\frac{\partial S}{\partial A} = 0$ . Logo:

$\begin{align*} \frac{\partial S}{\partial A} = &-2M^TX+A^TM^TM+M^TMA = 0\\ & 2M^TMA=2M^TX \end{align*}$

Por fim:

$A = (M^TM)^{-1}M^T X = M^\dagger X. \qquad \quad \qquad \square$

Obs. 1: Esse resultado é bem explicado no capítulo 3 do livro Algebra Linear and its Applications, de Gilbert Strang⁴[bibcite key=strang2006linear] e no capítulo 8 do livro Advanced Digital Signal Processing and Noise Reduction, de Saeed Vaseghi[bibcite key=vaseghi2008advanced]. Ficam como sugestões de leitura, para os mais interessados.

Obs. 2: Essa estruturação do problema é equivalente ao Método da Autocorrelação para a estimação dos coeficientes do preditor (sim, existem outros métodos para a estimação desses coeficientes). Porém a forma que utilizei aqui difere da estrutura que é mais comum de se ver nos artigos de predição linear e livros de processamento de sinais (um diferencial nesse aspecto é o livro do Vaseghi[bibcite key=vaseghi2008advanced]). Entretanto, trata-se de uma abordagem alternativa e mais intuitiva, sendo baseada em conceitos de álgebra linear.

Vamos testar essa estimação? (O sinal utilizado nesse trecho pode ser baixado nesse link: a.wav (361 downloads) )

Primeiro, vamos importar o arquivo que será utilizado:

[x,fa]=wavread('a.wav');

e em seguida, montaremos as matrizes X e M. A ordem $p=7$ foi utilizada aqui⁵.

p=7;
X=x';
M=zeros(length(x),p);
for n=1:length(x)
    if n<=p //Como valores menores que 1 são índices inválidos no scilab, é necessário fazer esse ajuste até a linha p 
        M(n,1:n-1)=x(n-1:-1:1);
    else
        M(n,:)=x(n-1:-1:n-p);

Por fim, podemos estimar $A$ :

A=inv(M'*M)*M'*X; //Pseudoinversa

e assim, podemos obter uma estimativa das amostras de $x$ ( $\hat x$ ) com base em suas amostras anteriores e, também, do erro de estimação $e$ . Na Figura 1 as primeiras 400 amostras dos sinais $x[n]$ , $\hat x[n]$ e $e[n]$ são ilustrados.

xhat=M*A;
plot(x(1:400))
plot(xhat(1:400),'r')
plot(e(1:400),'k')

**Figura 1** – Comparação entre amostras dos sinais $x[n]$ , $\hat x[n]$ e $e[n]$ , para o preditor de ordem 7.

Agora que já sabemos como estimar os coeficientes, vamos revisistar a equação a diferenças em (3):

$x[n]=\sum_{k=1}^p a_k x[n-k] \ + \ e[n]$

Ela parece familiar? Lembra alguma outra coisa estudada em processamento digital de sinais?

Se assumirmos que $e[n]$ seja a entrada de um filtro e $x[n]$ seja a sua saída, teremos aqui um filtro IIR. Dessa maneira, podemos considerar que o preditor linear vai nos permitir capturar as características de ressonância (polos) contidas no sinal analisado $x[n]$ . Então, podemos observar essas características de ressonância ao calcular a resposta em frequência do filtro do preditor. Para isso, aplicamos a Transformada Z em (3):

$X(z)=X(z)\left(\sum_{k=1}^p a_k z^{-k}\right)+E(z) ,$

$X(z)\left(1-\sum_{k=1}^p a_k z^{-k}\right) = E(z) ,$

de forma que a função de tranferência $H(z)=\frac{X(z)}{E(z)}$ é dada como segue:

(8) $\begin{equation*} H(z) = \frac{1}{\displaystyle 1-\sum_{k=1}^p a_k z^{-k}} . \end{equation*}$

Dessa maneira, podemos obter a resposta em frequência do filtro preditor (que vai representar as ressonâncias no sinal analisado). Devemos nos lembrar que para avaliar sobre o círculo unitário se assume que $z = e^{j\omega}$ , o que nos leva a representar (8) da seguinte maneira:

(9) $\begin{equation*} H(e^{j\omega}) = \frac{1}{\displaystyle 1-\sum_{k=1}^p a_k e^{-jk\omega}} . \end{equation*}$

Assim, para valores de $\omega$ entre $0$ e $\pi$ , calculamos o valor de $H(e^{j\omega})$ . Na Figura 2 se encontra representado o módulo de $H(e^{j\omega})$ em escala logaritmica.

w=0:0.01:%pi;
for i=1:length(w)
    H(i)=1/(1-exp(-%i*w(i)*[1:p])*A); // Cada posição de H representa o seu valor para cada valor de w
end
plot(w*fa/(2*%pi),log(abs(H)))
xgrid()

**Figura 2** – Resposta em frequência do preditor, para o sinal analisado, com ordem 7.

Antes de prosseguirmos, para maior comodidade, vamos encapsular os calculos do ajuste do preditor e de sua resposta em frequência em uma função no scilab.

function [A,Y,e,H,w]=predlinear(x,p,fa)
    X=x';
    M=zeros(length(x),p);
    for k=1:length(x)
        if k<=p
            M(k,1:k-1)=x([k-1:-1:1]);
        else
            M(k,:)=x([(k-1):-1:(k-p)]);
        end
    end
    A=inv(M'*M)*M'*X;
    Y=M*A;
    e=X-Y;
    w=0:0.01:%pi;
    for i=1:length(w)
        H(i)=1/(1-exp(-%i*w(i)*[1:p])*A);
    end
endfunction

Mas afinal, o que representa a ordem $p$ ? Qual é o efeito da alteração dessa ordem? Quando fazemos o ajuste dos coeficientes do preditor, minimizamos o erro quadrático utilizando o método dos mínimos quadrados. Mas o que será que acontece com o erro quadrático com a variação da ordem? Vamos testar? Os valores do erro quadrático para diferentes valores de $p$ são ilustrados na Figura 3.

for P=1:100 //Esse loop levar alguns segundos...
    [A,xc,e,H,w]=predlinear(x,P,fa);
    EQ(P)=sum(e^2);
end

plot(EQ)

Figura 3 – Erro quadrático para a ordem $p$ variando de 1 a 100.

Podemos observar que o erro quadrático segue diminuindo à medida em que a ordem aumenta. Então se o critério para a escolha da ordem fosse o menor erro quadrático, a maior ordem possível seria sempre escolhida.

Por outro lado, é válido levar em conta o efeito causado na resposta em frequência do preditor para a escolha da ordem. Se observarmos bem a equação da função de transferência em (8), podemos notar que o valor de $p$ define a quantidade polos que o preditor tenta estimar. Assim, quanto maior o número de polos, maior será a quantidade de ressonâncias observadas na resposta em frequência do preditor. Para observarmos isso, vamos plotar a resposta em frequência para $p = \{4,\ 8,\ 12,\ 16,\ 20\}$ . Os módulos das respostas em frequência, em dB ( $20\log_{10}\left(abs(H)\right)$ ), encontram-se ilustrados na Figura 4.

[A4,x4,e4,H4,w]=predlinear(x,4,fa);
[A8,x8,e8,H8,w]=predlinear(x,8,fa);
[A12,x12,e12,H12,w]=predlinear(x,12,fa);
[A16,x16,e16,H16,w]=predlinear(x,16,fa);
[A20,x20,e20,H20,w]=predlinear(x,20,fa);
[A24,x24,e24,H24,w]=predlinear(x,24,fa);

subplot(321)
plot(w*fa/(2*%pi),20*log10(abs(H4)),'linewidth',3)
xgrid()

subplot(322)
plot(w*fa/(2*%pi),20*log10(abs(H8)),'linewidth',3)
xgrid()

subplot(323)
plot(w*fa/(2*%pi),20*log10(abs(H12)),'linewidth',3)
xgrid()

subplot(324)
plot(w*fa/(2*%pi),20*log10(abs(H16)),'linewidth',3)
xgrid()

subplot(325)
plot(w*fa/(2*%pi),20*log10(abs(H20)),'linewidth',3)
xgrid()

subplot(326)
plot(w*fa/(2*%pi),20*log10(abs(H24)),'linewidth',3)
xgrid()

**Figura 4** – Efeito da ordem do preditor na resposta em frequência.

Podemos observar que, ao assumir valores muito baixos de $p$ , o preditor tende a gerar uma aproximação grosseira do espectro. Em contraste, valores altos da ordem nos impõem uma aproximação excessivamente detalhada do espectro. Assim, para a escolha da ordem $p$ é necessário se ter uma idéia prévia do sinal analisado e da quantidade de ressonâncias que esse deve conter.

Um exemplo interessante é a aplicação para sinais de voz (que por acaso é o que estamos testando). O sinal de voz é o resultado de uma excitação sonora, que ocorre na glote, ressonando em um tubo acústico que liga a glote até os lábios, denôminado trato vocal. As frequências de ressonância do trato vocal são denominadadas “formantes”. As 3 formantes mais relevantes oscilam em torno de 500, 1500 e 2500 Hz (variando de acordo com a vogal e com características anatômicas e fisiólogica de quem fala), respectivamente. Logo, uma ordem adequada para o preditor deve garantir uma estimação de uma resposta em frequência que englobe, ao menos, picos em valores razoáveis para as três primeiras formantes.

Uma outra forma interessante para observar a adequação espectral do preditor, para uma ordem $p$ qualquer, é a comparação da resposta em frequência estimada com o espectro da transformada discreta de Fourier. Uma comparação entre o espectro da DFT e a resposta em frequência para $p = \{4,\ 16,\ 64\}$ é ilustrada na Figura 5. Nela é possível ver a influência da ordem no ajuste da resposta em frequência ao espectro de Fourier.

[A4,x4,e4,H4,w]=predlinear(x,4,fa);
[A16,x16,e16,H16,w]=predlinear(x,16,fa);
[A64,x64,e64,H64,w]=predlinear(x,64,fa);

Fx=fft(x); // Função do scilab que executa a FFT -> algoritmo que realiza de maneira rápida a DFT

plot(linspace(0,fa/2,ceil(length(Fx)/2)),20*log10(abs(Fx(1:ceil(length(Fx)/2)))),'m')
plot(w*fa/(2*%pi),20*log10(abs(H4)),'linewidth',5)
plot(w*fa/(2*%pi),20*log10(abs(H16)),'g','linewidth',5)
plot(w*fa/(2*%pi),20*log10(abs(H64)),'k','linewidth',5)
xgrid()

**Figura 5** – Comparação entre espectro de Fourier e as respostas em frequências de preditores de diferentes ordens.

Outras considerações válidas para escolher a ordem do preditor podem ser encontradas na literatura. Para os interessados em se aprofundar mais no assunto, boas opções de leitura, além do capítulo 8 do Vaseghi [bibcite key=vaseghi2008advanced], são o capítulo 8 do livro Digital Signal Processing of Speech Signals, de Rabiner e Schafer[bibcite key=Rabiner1978], e os artigos de John Makhoul Linear Prediction: A Tutorial Review[bibcite key=makhoul1975linear] e Spectral Analysis of Speech by Linear Prediction[bibcite key=makhoul1973spectral].

Uma vez decidida a ordem e com os coeficientes estimados, temos um filtro com um conjunto de ressonâncias que representam o trato vocal do indíviduo que pronunciou a vogal. Que tal se usássemos esse filtro então? Se colocássemos, por exemplo, um ruído branco na entrada desse filtro, deveríamos esperar que o sinal resultante adquira as características ressonantes da vogal emitida no sinal original. Ou seja, o som do sinal resultante deverá se parecer com a vogal emitida no sinal. Vamos testar?

Utilizando o sinal que vinhamos utilizando (equivalente à vogal /a/), tudo que precisaremos fazer é substituir o sinal $e[n]$ pelo sinal desejado em (3). Assim, utilizamos os coeficientes obtidos pela estimação do preditor para filtrar o sinal aplicado à entrada. Para facilitar esse processo, criei uma pequena função para implementar o filtro IIR com base apenas nos coeficientes $a_k$ .

function y=iirFiltro(x,a)
    y=zeros(1,length(x));
    memo=zeros(1,length(a));
    for n=1:length(x)
        y(n)=x(n)+memo*a;
        memo=[y(n),memo];
        memo($)=[];
    end
    y=y/max(abs(y));
endfunction

Vamos testar essa implementação para um ruído branco e para um trem de impulsos como entradas dos filtros. Aqui utilizaremos um preditor de ordem $p=16$ .

[A,xc,e,H,w]=predlinear(x,16,fa);

Para um ruído branco:

r=rand(length(x),1);
xr=iirFiltro(r,A);

sound(r,fa)

sound(xr,fa)

sound(x,fa)

Na Figura 6 o efeito da filtragem no espectro de Fourier é ilustrado.

Fr=fft(r);
Fxr=fft(xr);
plot(linspace(0,fa/2,ceil(length(Fr)/2)),20*log10(abs(Fr(1:ceil(length(Fr)/2)))),'k')
plot(linspace(0,fa/2,ceil(length(Fxr)/2)),20*log10(abs(Fxr(1:ceil(length(Fxr)/2)))),'m')
plot(w*fa/(2*%pi),20*log10(abs(H)),'linewidth',5)

**Figura 6** – Efeito espectral da filtragem para o ruído branco.

Para um trem de Impulsos com frequência fundamental de 200Hz:

I=zeros(1,length(x));
f0=200;
I(1:fa/f0:length(I))=1;
xi=iirFiltro(I,A);

sound(r,fa)

sound(xr,fa)

sound(x,fa)

Na Figura 7 o efeito da filtragem no espectro de Fourier é ilustrado.

Fxi=fft(xi);
Fi=fft(I);
plot(linspace(0,fa/2,ceil(length(Fi)/2)),20*log10(abs(Fi(1:ceil(length(Fi)/2)))),'k','linewidth',2)
plot(linspace(0,fa/2,ceil(length(Fxi)/2)),20*log10(abs(Fxi(1:ceil(length(Fxi)/2)))),'m','linewidth',2)
plot(w*fa/(2*%pi),20*log10(abs(H)),'linewidth',5)

**Figura 7** – Efeito espectral da filtragem para o trem de impulsos.

Bem legal, né?

Por fim, seguem algumas considerações para finalizarmos:

– Se analisarmos (3) com “outros olhos” e assumirmos $e[n]$ como saída e $x[n]$ como entrada, teremos um filtro FIR ( $e[n] = x[n] - \sum_{k=1}^p x[n-k]$ ). Ora, esse filtro terá a resposta em frequência inversa ao filtro IIR (i.e., $G(z)={H(z)}^-1$ , onde $G(z)$ é a resṕosta em frequência do filtro FIR). Assim, o sinal $e[n]$ será o resultado da remoção do efeito das ressonâncias estimadas do sinal original (troca de polos por zeros nos mesmos locais). Esse processo é denominado filtragem inversa.

– Ao longo do desenvolvimento desse texto, para simplificar a análise, o sinal foi sempre utilizado por inteiro. Entretanto, isso tende a prejudicar os resultados, pois o trato vocal tende se alterar com o tempo (não estamos lidando com um sistema invariante no tempo). Uma boa prática para melhores resultados é o uso de segmentos menores do sinal, janelados no tempo. Normalmente, para sinais de voz, uma janela com 20 milisegundos de duração é recomendada para a análise, pois nesse tempo não ocorrem alterações significativas no trato vocal.

– O capítulo 8 do Vaseghi[bibcite key=vaseghi2008advanced] é uma boa leitura complementar, que apresenta os conceitos aqui apresentados (e vai além) de uma maneira bem intuitiva. As demais sugestões de leitura aqui apresentadas, apesar de não tão intuitivos quanto o Vaseghi, são boas opções para quem tem o interesse de se aprofundar no tema.

Em caso de dúvidas, sinta-se a vontade para comentar (no final da página) e assim que possível, tentarei esclarecer.

Referências / Sugestões de Leitura

[/bibshow]

Sobre o Autor

Jônatas Cruz Santos

jonatascrs@live.com | CV Lattes

Mestrando em Engenharia Elétrica e Bacharel em Engenharia Eletrônica pela Universidade Federal de Sergipe. Interessado nas áreas de Modelagem Computacional, Processamento Digital de Sinais, Processamento Digital da Fala e Reconhecimento de Padrões. Atualmente desenvolve pesquisa no estudo de modelos computacionais dos mecanismos de produção vocal.

Áreas de Pesquisa: Modelagem Computacional, Processamento Digital de Sinais, Reconhecimento de Padrões, Processamento de Sinais de Voz e Fala, Qualidade Vocal.

Optou-se por utilizar os índices começando em 1 para que haja compatibilidade com a notação do Scilab
O uso de $\hat x$ representa uma estimativa de $x$
$M^\dagger = \text{Pseudoinversa de M}$
Não confundir Dr. Strang com Dr. Strange
Por sugestão completamente aleatória de Molininha

Biochaves

Introdução à Predição Linear (com aplicação em sinais de voz)

– Chamamos um vidente?

– Usamos o Olho de Agamotto?

– Sentamos e choramos?

Referências / Sugestões de Leitura

Jônatas Cruz Santos

1 thought on “Introdução à Predição Linear (com aplicação em sinais de voz)”

Deixe uma resposta Cancelar resposta

Posts recentes

Categorias

Se inscreva na nossa Newsletter!