Ao final desta unidade, o aluno
Antes de começar a discussão desta unidade, o aluno deve estar tranquilo com os seguintes conceitos:
Se um sinal $y[n]$ é obtido pela convolução entre uma entrada $x[n]$ e uma resposta ao impulso $h[n]$, então temos que sua Transformada de Fourier é calculada pela multiplicação entre as Transformadas de Fourier da entrada e do impulso, ou seja: $$ y[m] = \sum _{m=-M}^M x[m]h[n-m] \iff Y[k] = X[k]H[k].$$
Para calcular uma somatória num programa de computador, é possível utilizar um laço for. Assim, o resultado da expressão: $$y[n] = \sum _{m=0}^n \frac{x[m]}{n}$$
pode ser calculado pelo programa:
In [2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
n_max = 100 # Este valor é arbitrário
y = np.zeros(n_max)
x = np.random.random(n_max)
for n in xrange(n_max):
for m in xrange(n):
y[n] += x[m]/float(n)
plt.figure()
plt.plot(x)
plt.plot(y)
plt.show()
Cordas e tubos vibrantes oscilam em diversas frequências, todas múltiplas de uma frequência fundamental $f_0$. Por este motivo, o sinal de áudio resultante desta vibração pode ser razoavelmente bem descrito por uma Série de Fourier, na forma:
$$x(t) = \sum_{m=1}^M a_m \cos(2 \pi m f_0 t + \phi_m).$$Um sinal $x(t)$ é periódico com período $\tau$ se $x(t) = x(t-\tau)$ para todo $t$.
A auto-correlação de um sinal é um sinal de tamanho igual ao sinal original, calculado por: $$r_{x}[k] = \sum_{n=-\infty}^{\infty} x[n]x[n+k].$$
A auto-correlação representa a energia de um sinal multiplicado por ele mesmo e deslocado $k$ amostras. Trata-se de um sinal que é máximo para $k=0$ e que decai com o aumento de $k$.
Na prática, sinais $x[n]$ sempre têm uma duração conhecida. Por isso, os limites da somatória sempre são modificados de forma a contemplar essa duração conhecida.
Vejamos o que acontece quando calculamos algumas auto-correlações:
In [3]:
def acc(x):
"""Calcula a auto-correlacao de um sinal x[n]"""
r = np.zeros(x.shape)
for k in xrange(len(r)):
for n in xrange(len(x)):
if (n+k) < len(x): # Esta condicional evita acessar elementos x[m] para m < 0
r[k] += x[n] * x[n+k]
return r
rx = acc(x)
ry = acc(y)
plt.figure()
plt.plot(rx)
plt.plot(ry)
plt.show()
Veja que, neste caso, há progressivamente menos amostras não-nulas no sinal de auto-correlação. Poderíamos corrigir isso utilizando uma janela deslizante com a metade do tamanho de $x[n]$, e restringindo nosso cálculo à metade do comprimento de $x[n]$.
In [4]:
def acc2(x):
"""Calcula a auto-correlacao de um sinal x[n]"""
r = np.zeros(len(x)/2)
for k in xrange(len(r)):
for n in xrange(len(r)):
if (n+k) > 0: # Esta condicional evita acessar elementos x[m] para m < 0
r[k] += x[n] * x[n+k]
return r
rx = acc2(x)
ry = acc2(y)
plt.figure()
plt.plot(rx)
plt.plot(ry)
plt.show()
Uma propriedade importante da auto-correlação é que se $x(t)$ é periódico então sua autocorrelação $r_x(t)$ também é periódica, com mesmo período. Veja o que acontece quando geramos uma onda quadrada e calculamos sua auto-correlação:
In [5]:
t = np.linspace(0, 1, 100)
x = np.sign( np.cos(2*np.pi*10*t))
r = acc2(x)
plt.figure()
plt.plot(t,x)
plt.xlabel('Tempo (s)')
plt.show()
plt.figure()
plt.plot(t[0:len(r)], r)
plt.xlabel('Deslocamento (s)')
plt.show()
Verificamos que $r_x$ apresenta picos espaçados de 0.1s, o que corresponde ao período fundamental de nosso sinal $x[n]$. Portanto, precisamos de algum algoritmo para encontrar o primeiro pico da autocorrelação que não seja o correspondente ao deslocamento zero. Embora existam muitas propostas na literatura, é possível fazer o seguinte:
In [6]:
def upsample_subtract(r):
r2 = np.zeros(len(r)/2)
for n in xrange(len(r2)):
r2[n] = r[n]-r[n/2]
return r2
r2 = upsample_subtract(r)
plt.figure()
plt.plot(t[0:len(r2)], r2)
plt.xlabel('Deslocamento (s)')
plt.show()
r2 *= np.linspace(1, 0, len(r2))
t0 = np.argmax(r2)
plt.figure()
plt.plot(t[0:len(r2)], r2)
plt.plot(t[t0], r2[t0], 'ro')
plt.xlabel('Deslocamento (s)')
plt.show()
Neste ponto, portanto, temos um algoritmo capaz de estimar o período fundamental em um pequeno trecho contendo uma janela de um sinal periódico. Vamos melhorá-lo.
Para o cálculo de cada amostra da auto-correlação, é preciso realizar $N$ somas, onde $N$ é o tamanho do sinal original $x[n]$. Como são $N$ amostras, então calcular uma auto-correlação precisa de $O(N^2)$ cálculos, o que pode se tornar um problema. Porém, podemos re-escrever a auto-correlação como: $$r_{x}[k] = \sum_{n=-\infty}^{\infty} x[n]x[n+k]. = \sum_{n=-\infty}^{\infty} x[n]x[n-(-k)],$$ mostrando que o cálculo da auto-correlação, na verdade, é uma convolução de um sinal com ele mesmo. Portanto, a auto-correlação pode ser calculada no domínio da frequência, reduzindo a sua complexidade computacional para $O(n \log_2(n))$, dando origem ao seguinte algoritmo para calculo do periodo fundamental:
In [7]:
def acc_fd(x):
X = np.abs(np.fft.fft(x))
r = np.real(np.fft.ifft(X*X))
return r
r = acc_fd(x)
plt.figure()
plt.plot(t,x)
plt.xlabel('Tempo (s)')
plt.show()
plt.figure()
plt.plot(t[0:len(r)], r)
plt.xlabel('Deslocamento (s)')
plt.show()
r2 = upsample_subtract(r)
plt.figure()
plt.plot(t[0:len(r2)], r2)
plt.xlabel('Deslocamento (s)')
plt.show()
r2 *= np.linspace(1, 0, len(r2))
t0 = np.argmax(r2)
plt.figure()
plt.plot(t[0:len(r2)], r2)
plt.plot(t[t0], r2[t0], 'ro')
plt.xlabel('Deslocamento (s)')
plt.show()
Assim, podemos consolidar nosso algoritmo numa única função, na forma:
In [8]:
def t0_acc(x):
"""Retorna o periodo fundamental de x, em amostras"""
X = np.abs(np.fft.fft(x))
r = np.real(np.fft.ifft(X*X))
r2 = np.zeros(len(r)/2)
for n in xrange(len(r2)):
r2[n] = r[n]-r[n/2]
r2 *= np.linspace(1, 0, len(r2))
t0 = np.argmax(r2)
return t0
Até o momento, aplicamos o detector de frequências fundamentais apenas em sinais fictícios. Gostaríamos de aplicá-lo em sinais de áudio reais. Para fazer isso, é preciso lembrar que uma mesma frequência fundamental não será adequada para descrever um sinal de áudio inteiro, mas, dentro de uma pequena janela (de 23ms, por exemplo), é pouco provável que a frequência fundamental mude. Também, precisamos lembrar que o algoritmo usando autocorrelação requer uma janela que seja pelo menos duas vezes maior que o período fundamental que será medido. Assim, temos:
In [12]:
%matplotlib inline
import numpy as np
import scipy.io.wavfile
import matplotlib.pyplot as plt
import IPython.lib.display as display
fname = 'audio/piano.wav'
rate, data = scipy.io.wavfile.read(fname)
data = data.astype(np.float)
t0 = 0
window = 1024
step = 512
f0 = []
while t0+window < len(data):
per_fun = t0_acc(data[t0:t0+window])
if per_fun > 2:
per_fun /= float(rate)
f0.append(1./per_fun)
else:
f0.append(0)
t0 += step
t = np.linspace(0, len(data)/rate, len(f0))
plt.plot(t,f0)
plt.ylabel('Frequencia Fundamental (Hz)')
plt.xlabel('Tempo (s)')
plt.show()
display.Audio(fname)
Out[12]:
Podemos perceber que:
Vamos ver o que acontece ao processarmos um glissando:
In [10]:
fname = 'audio/chirp.wav'
rate, data = scipy.io.wavfile.read(fname)
data = data.astype(np.float)
t0 = 0
window = 1024
step = 512
f0 = []
while t0+window < len(data):
per_fun = t0_acc(data[t0:t0+window])
if per_fun > 2:
per_fun /= float(rate)
f0.append(1./per_fun)
else:
f0.append(0)
t0 += step
t = np.linspace(0, len(data)/rate, len(f0))
plt.plot(t,f0)
plt.ylabel('Frequencia Fundamental (Hz)')
plt.xlabel('Tempo (s)')
plt.show()
display.Audio(fname)
Out[10]:
O sistema parece ter se comportado um pouco melhor. Porém, veja que, ao invés de estimar uma curva crescente, na verdade estimamos uma "escada". Isso acontece porque nosso sistema de detecção tem precisão finita.
Por fim, vamos tentar rodar o algoritmo em um solo de Blues. Adicionalmente, para esta análise, converteremos nossas frequências para números MIDI:
In [11]:
fname = 'audio/bbking.wav'
rate, data = scipy.io.wavfile.read(fname)
data = data.astype(np.float)
t0 = 0
window = 1024
step = 512
f0 = []
while t0+window < len(data):
per_fun = t0_acc(data[t0:t0+window])
if per_fun > 2:
per_fun /= float(rate)
f0.append(69+12*np.log2(1./(per_fun*440)))
else:
f0.append(0)
t0 += step
t = np.linspace(0, len(data)/rate, len(f0))
plt.plot(t,f0)
plt.ylabel('Numero MIDI')
plt.xlabel('Tempo (s)')
plt.show()
display.Audio(fname)
Out[11]:
Neste caso, exceto pelos comportamentos anômalos nos ataques (que são regiões em que tipicamente não vemos comportamento harmônico), o algoritmo parece ter retornado informações bastante precisas.
In [ ]: