[Could not find the bibliography file(s)

Suppose that, by sampling a random signal, we have obtained the following sequence:

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

and, using psychic (or mathematical, preferably) powers, we would like to estimate its next sample. How should we proceed?

– Should we call a Fortune Teller?
– Should we use Eye of Agamotto?
– Should we sit and cry?

Well, some digital signal processing concepts and a good shot of linear algebra might give us more viable alternatives to solve this problem. First of all, to use a standard and easy notation, let's call this sequence x[n], such that x[1]=2.0, x[2]=3.5, etc.1 An interesting idea would be to use a linear combination of the previous samples. Thus, we could estimate the next sample of this signal by using, for example, 3 previous samples as follows:

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

Essas constantes (a_1a_2a_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 maneira2:

    \[ \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_1a_2a_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 maneira3:

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

ou

(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 Strang4[?] e no capítulo 8 do livro Advanced Digital Signal Processing and Noise Reduction, de Saeed Vaseghi[?]. 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[?]). 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 (91 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 aqui5.

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 [?], são o capítulo 8 do livro Digital Signal Processing of Speech Signals, de Rabiner e Schafer[?], e os artigos de John Makhoul Linear Prediction: A Tutorial Review[?] e Spectral Analysis of Speech by Linear Prediction[?].

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[?] é 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

About the Author

Graduate Student in Electrical Engineering (Masters Degree) and Graduated in Electronic Engineering at Federal University of Sergipe - Brazil. Interested in Computational Modelling, Digital Signal Processing, Digital Processing of Speech and Patterns Recognition. Nowadays develops a research, computational models of voice production mechanisms.

Research Topics: Computational Modelling, Digital Signal Processing, Pattern Recognition, Voice and Speech Signal Processing, Voice Quality.

  1. It was chosen to use index beginning at 1, for SciLab compatibility.
  2. O uso de \hat x representa uma estimativa de x
  3. M^\dagger = \text{Pseudoinversa de M}
  4. Não confundir Dr. Strang com Dr. Strange 
  5. Por sugestão completamente aleatória de Molininha

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

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.

en_GBEnglish (UK)
pt_BRPortuguês do Brasil en_GBEnglish (UK)