Exemplos em Octave (síntese e análise) da aula de 05/08/14

#!/usr/bin/octave
#
# Este é um código escrito para a linguagem Octave (www.octave.org)
#
# Funciona para a versão 3.0.1 mas deve funcionar em outras versões próximas.
#
# Este arquivo pode ser executado (em Linux) no interpretador de comandos,
# ou dentro do interpretador do octave (digitando aula1exemplos), ou ainda
# copiando e colando cada porção para o interpretador do octave.
#
# Duas observações importantes:
#
# 1) Para compreender o modo como os exemplos foram construídos, apenas
#    rodar o arquivo inteiro no octave não é suficiente, é preciso
#    olhar o código com cuidado.
#
# 2) Nos exemplos abaixo não são fornecidas muitas explicações de porque
#    as experiências têm estes resultados e não outros, e como se pode
#    prever ou controlar estes resultados. O ÚNICO OBJETIVO DESTA AULA
#    É ESTIMULAR A CURIOSIDADE SOBRE ESTES PROCESSOS, e se ela conseguir
#    o que pretende, terá sido imensamente produtiva.
#
# Os exemplos utilizam o comando aplay da biblioteca ALSA (Linux)
# para tocar arquivos wav. Em outras plataformas, você deve substituir
# este comando por qualquer aplicativo de linha de comando que toque
# arquivos de áudio.
#

# Faz o terminal do Octave imprimir as mensagens imediatamente
page_screen_output(false);

disp("#");
disp("# Este primeiro exemplo cria um sinal senoidal (que");
disp("# corresponde a um movimento harmônico simples) com");
disp("# frequência de oscilação de 440Hz, usando uma");
disp("# frequência de amostragem de 44100 valores por");
disp("# segundo, durante 2 segundos, com amplitude de");
disp("# oscilação máxima (=1 por convenção).");
disp("#");
disp("Pressione qualquer tecla...");
input("");

# cria um sinal senoidal e armazena em um vetor (coluna)
sinal1 = sinetone(440,44100,2,1);
disp("# Eis os primeiros 10 valores deste vetor:");
disp(sinal1(1:10));
disp("Pressione qualquer tecla...");
input("");
disp("# Constrói o gráfico com os primeiros 1000 valores,");
disp("# utilizando gnuplot");
plot(sinal1(1:1000));
disp("#");
disp("Pressione qualquer tecla...");
input("");
disp("# Escreve o sinal em um arquivo de áudio");
disp("# e toca o arquivo usando aplay");
wavwrite(sinal1,44100,"temp.wav");
system("aplay temp.wav");
disp("#");
disp("Pressione qualquer tecla...");
input("");

#
# Você pode obter a documentação dos comandos
# sinetone, disp, plot, wavwrite e system digitando
# doc nomedocomando






# Uma série harmônica e um exemplo de aliasing.

disp("# Neste exemplo construiremos um sinal que varre");
disp("# os harmônicos de 440Hz (440, 880, 1320, ...)");
disp("# para ilustrar a diferença que pode fazer");
disp("# a escolha da taxa de amostragem na");
disp("# representação digital.");
disp("#");
disp("Pressione qualquer tecla...");
input("");

disp("# Toca a série com SR=44100Hz");

SR = 44100;

# cria um vetor vazio
x=[];
# e vai anexando cada um dos 16 primeiros harmônicos
for k=1:16
  x = [x;sinetone(440*k, SR, 0.5, 1)];
endfor
# escreve no arquivo de áudio e toca
wavwrite(x,SR,"temp.wav");
system("aplay temp.wav");
disp("#");
disp("Pressione qualquer tecla...");
input("");

disp("# Toca a série com SR=8000Hz");

SR = 8000;

x=[];
for k=1:16
  x = [x;sinetone(440*k, SR, 0.5, 1)];
endfor
wavwrite(x,SR,"temp.wav");
system("aplay temp.wav");

disp("#");
disp("VOLTAR PARA A APRESENTAÇÃO!");
disp("Pressione qualquer tecla...");
input("");





# define a taxa de amostragem (Sample Rate) como variável global
SR=44100;


# a função abaixo serve para poupar digitação, chamando para
# cada sinal o comando plot, gravando o sinal em arquivo wav
# e reproduzindo o arquivo.
function plotnplay(SR,s)
  # chama gnuplot. O comando drawnow() força a atualização da janela,
  # que normalmente só seria desenhada no retorno ao prompt do octave
  plot(s(1:1000));drawnow();
  # escreve arquivo wav temporário
  wavwrite(s,SR,"temp.wav");
  # toca arquivo usando ALSA
  system("aplay temp.wav");
  # apaga arquivo temporário
  system("rm temp.wav");
endfunction






disp("# ");
disp("# O próximo exemplo cria uma mistura (soma) de 3");
disp("# sinais senoidais de frequências diferentes, que");
disp("# formam um acorde maior A4 - C#5 - E5 (estes são");
disp("# harmônicos da nota A2 que corresponde à frequência");
disp("# de 110Hz)");
disp("#");
disp("Pressione qualquer tecla...");
input("");

# cria a mistura das 3 senóides. Observe que as amplitudes foram
# ajustadas para 0.3 a fim de impedir que o sinal "estourasse"
# o limite da faixa [-1...+1]
sinal2 = sinetone(440,SR,2,0.3)+sinetone(550,SR,2,0.3)+sinetone(660,SR,2,0.3);
plotnplay(SR,sinal2);


disp("# ");
disp("# O exemplo a seguir cria outra mistura de 3");
disp("# senóides, porém com frequências aleatórias. São");
disp("# sorteados 3 valores entre 220Hz e 440Hz.");
disp("#");
disp("Pressione qualquer tecla...");
input("");

# sorteia vetor coluna com 3 frequências entre 220 e 440
freq = round(220*rand(3,1))+220;
# cria a mistura dos 3 sinais senoidais
sinal3 = sinetone(freq(1),SR,2,0.3)+sinetone(freq(2),SR,2,0.3)+sinetone(freq(3),SR,2,0.3);
plotnplay(SR,sinal3);

disp("# ");
disp("# O objetivo deste exemplo na realidade é mostrar");
disp("# como podemos descobrir as frequências que compõem");
disp("# o sinal utilizando um truque trigonométrico");
disp("# simples, chamado de ortogonalidade dos senos e");
disp("# cossenos. Esta propriedade é consequência da");
disp("# expressão cos(a)cos(b) = 0.5*[cos(a+b)+cos(a-b)],");
disp("# que pode ser reformulada para vetores de amostras");
disp("# com frequências f e g como cos(f*n)cos(g*n) =");
disp("# 0.5*[cos((f+g)*n)+cos((f-g)*n)].  Em outras");
disp("# palavras, o produto dos sinais senoidais com");
disp("# frequências f e g é idêntico à média de outros");
disp("# dois sinais senoidais, com frequências f+g e");
disp("# f-g. Ao somar todas as amostras do produto dos");
disp("# sinais originais, teremos o mesmo resultado que ao");
disp("# somar as amostras dos sinais senoidais com");
disp("# frequências f+g e f-g.");
disp("# ");
disp("Pressione qualquer tecla...");
input("");
disp("# ");
disp("# Há dois casos possíveis:");
disp("# ");
disp("# 1) f e g são diferentes: neste caso os valores");
disp("# positivos das senoides f+g e f-g serão cancelados");
disp("# pelos valores negativos das mesmas senoides, e a");
disp("# soma será ZERO");
disp("# ");
disp("# 2) f e g são iguais: neste caso a expressão");
disp("# cos((f-g)*n) é sempre 1 (pois f-g=0), e a soma");
disp("# destes valores será INFINITO (a outra senoide, de");
disp("# frequência f+g terá o mesmo cancelamento de");
disp("# valores do caso 1).");
disp("# ");
disp("Pressione qualquer tecla...");
input("");
disp("# ");
disp("# Por exemplo, as senóides com frequências 440 e");
disp("# 110Hz são ortogonais, sendo a soma acima quase");
disp("# ZERO:");
# calcula o produto ponto-a-ponto dos vetores de amostras,
# e soma todas as posições do vetor resultante
printf("sum(sinetone(440,SR,2,1).*sinetone(100,SR,2,1))=%g\n",sum(sinetone(440,SR,2,1).*sinetone(100,SR,2,1)));
disp("# ");
disp("# (a razão da soma não ser exatamente ZERO é que");
disp("# estamos trabalhando com senóides truncadas num");
disp("# vetor de tamanho finito)");
disp("# ");
disp("# Um exemplo análogo com as frequências próximas");
disp("# 440 e 441Hz:");
printf("sum(sinetone(440,SR,2,1).*sinetone(441,SR,2,1))=%g\n",sum(sinetone(440,SR,2,1).*sinetone(441,SR,2,1)));
disp("# ");
disp("Pressione qualquer tecla...");
input("");
disp("# ");
disp("# Para ilustrar o caso 2, tomamos o produto de duas");
disp("# senóides com frequência 440Hz. Todos os valores do");
disp("# vetor produto são positivos, e a soma é muito");
disp("# grande (não é INFINITO pois o vetor tem tamanho");
disp("# finito):");
disp("#");
printf("sum(sinetone(440,SR,2,1).*sinetone(440,SR,2,1))=%g\n",sum(sinetone(440,SR,2,1).*sinetone(440,SR,2,1)));
disp("Pressione qualquer tecla...");
input("");


disp("# ");
disp("# Voltando ao caso das 3 senóides misteriosas,");
disp("# podemos procurar suas frequências por tentativa e");
disp("# erro, ou podemos fazer uma busca exaustiva num");
disp("# conjunto razoavelmente grande de frequências");
disp("# possíveis. Vamos ilustrar esta busca considerando");
disp("# todas as frequências inteiras entre 1 e 1000Hz,");
disp("# para cada uma delas calcular o produto da senóide");
disp("# correspondente com o sinal anterior e fazer a soma");
disp("# dos elementos do vetor. Para cada frequência");
disp("# teremos um valor da soma, e estes valores serão");
disp("# colocados em um gráfico: (aguarde um instante!)");
disp("#");

# gera um vetor coluna de tamanho 1000, inicializado com zeros
soma = zeros(1000,1);
# para cada frequência f entre 1 e 1000, calcula o produto com a
# senóide correspondente, e guarda a soma dos valores em soma(f)
for f=1:1000
  soma(f)=sum(sinetone(f,SR,2,1).*sinal3);
endfor
# mostra o gráfico
plot(soma);
# Para identificar as posições do vetor onde a soma foi grande,
# criamos um vetor de 0's e 1's para marcar estas frequências
# (através da expressão (soma>1)), depois multiplicamos ponto-a-ponto
# pelo vetor de índices de 1 a 1000 (1:1000) e recortamos apenas
# os valores diferentes de zero:
disp("# As frequências em destaque no gráfico são:");
disp("#");
disp(nonzeros((soma>1).*(1:1000)'));
disp("# As frequências anteriormente sorteadas eram:");
disp("#");
disp(freq);

disp("Pressione qualquer tecla...");
input("");


disp("# ");
disp("# Para mostrar a relação do que acabamos de fazer");
disp("# com o que se denomina normalmente de Análise");
disp("# Harmônica (ou de Fourier), veremos a seguir o");
disp("# gráfico do espectro de magnitude do sinal");
disp("# anterior, calculado usando a função fft() do");
disp("# octave. A função abs() calcula a magnitude ou");
disp("# valor absoluto (pois o resultado da fft() é");
disp("# complexo), e o gráfico mostra a porção inicial");
disp("# deste espectro. (os índices não têm os mesmos");
disp("# valores das frequências do exemplo anterior, mas");
disp("# correspondem a eles através de uma expressão");
disp("# simples).");
plot(abs(fft(sinal3))(1:2000));
disp("VOLTAR PARA A APRESENTAÇÃO!");
disp("Pressione qualquer tecla...");
input("");


disp("# ");
disp("# A seguir um fragmento de áudio é obtido de um");
disp("# arquivo wav e armazenado em um vetor do octave,");
disp("# para futuras manipulações:");
disp("#");
disp("Pressione qualquer tecla...");
input("");

sinal4 = wavread("exemplo.wav");
sinal4 = sinal4(1:5*44100,1);
plotnplay(SR,sinal4);





disp("# ");
disp("# Por manipulações, entendemos quaisquer operações");
disp("# feitas sobre o arquivo cujo resultado seja também");
disp("# um vetor de amostras com valores na faixa");
disp("# [-1...+1]. Por exemplo, podemos tomar médias de");
disp("# amostras adjacentes, definindo a partir de um");
disp("# sinal x um novo sinal y=0.5*(x+x(n+1)):");
disp("#");
disp("Pressione qualquer tecla...");
input("");

# obtém o comprimento do vetor (número de amostras total)
n=length(sinal4);
# calcula a fórmula acima "em bloco", ou seja, somando os
# subvetores do sinal 4 com índices de 1 a n-1 (bloco 1)
# e de 2 a n (bloco 2), e dividindo por 2
x = (sinal4(1:n-1)+sinal4(2:n))/2;
plotnplay(SR,x);



disp("# ");
disp("# Alguma coisa mudou no som, mas provavelmente a");
disp("# mudança é muito sutil para ser percebida. Podemos");
disp("# tornar esta mudança mais evidente fazendo médias");
disp("# não apenas de duas amostras vizinhas, mas de");
disp("# conjuntos de k amostras, através da fórmula");
disp("# y=(x+x(n+1)+x(n+2)+...+x(n+k-1))/k.");
disp("#");
disp("Pressione qualquer tecla...");
input("");

# a função abaixo devolve um sinal obtido de s substituindo
# cada amostra por uma média de k amostras consecutivas.
# A resposta terá tamanho n-k+1, onde n é o tamanho da entrada
function m = media(s,k)
  # n é o tamanho da entrada
  n = length(s);
  # m (a saída) é inicializado com zeros
  m = zeros(n-k+1,1);
  # a média é calculada em bloco outra vez, desta vez somando-s
  # os k blocos de índices (1:n-k+1), (2:n-k+2), ..., (k:n)
  for i=1:k
    m += s(i:n-k+i)/k;
  endfor
endfunction


disp("# ");
disp("# Para k=2, temos:");
plotnplay(SR,media(sinal4,2));
disp("# ");
disp("# Para k=10, temos:");
plotnplay(SR,media(sinal4,10));
disp("# ");
disp("# e para k=50, temos:");
disp("#");
plotnplay(SR,media(sinal4,50));

disp("# ");
disp("# A mudança na qualidade do som deve ser");
disp("# perceptível, mas a verdadeira questão é: porquê");
disp("# isso aconteceu, como poderíamos prever e como");
disp("# poderíamos controlar este resultado. É disso que");
disp("# trata a teoria de filtros para processamento de");
disp("# sinais.");
disp("#");
disp("Pressione qualquer tecla...");
input("");





disp("# ");
disp("# Outro exemplo simples de manipulação corresponde a");
disp("# calcular diferenças entre amostras (seguidas ou");
disp("# não). Estas diferenças correspondem a versões");
disp("# discretizadas da derivada do sinal, como na");
disp("# expressão y = x(n+1)-x, ou ainda y =");
disp("# x(n+k)-x. Para garantir que o resultado");
disp("# permanece na faixa [-1...+1], as expressões devem");
disp("# ser multiplicadas por 0.5.");
disp("#");
disp("Pressione qualquer tecla...");
input("");

# a função abaixo calcula as diferenças de ordem k do sinal de entrada
function d = diferenca(s,k)
  # n é o tamanho da entrada
  n = length(s);
  # a operação mais uma vez é feita em bloco
  d = s(k+1:n)-s(1:n-k);
endfunction

disp("# ");
disp("# Considerando a fórmula y = 0.5*(x(n+1)-x)");
disp("# temos este resultado.");
plotnplay(SR,diferenca(sinal4,1));
disp("# ");
disp("# Considerando a fórmula y = 0.5*(x(n+2)-x)...");
plotnplay(SR,diferenca(sinal4,2));
disp("# ");
disp("# Considerando a fórmula y = 0.5*(x(n+10)-x)...");
plotnplay(SR,diferenca(sinal4,10));
disp("# ");
disp("# e finalmente a fórmula y = 0.5*(x(n+50)-x)...");
plotnplay(SR,diferenca(sinal4,50));
disp("#");





disp("# Provavelmente as mudanças são todas bem");
disp("# perceptíveis, mas não seguem uma lógica parecida");
disp("# com a do exemplo anterior, no sentido que aumentar");
disp("# o parâmetro k não parece exatamente intensificar o");
disp("# efeito que foi percebido para k=1...");
disp("# ");
disp("Pressione qualquer tecla...");
input("");





disp("# ");
disp("# A sequência de manipulações abaixo considera");
disp("# produtos de sinais com valores na faixa [-1...+1]");
disp("# (cujo resultado está sempre na faixa [-1...+1]).");
disp("# ");
disp("# Serão retomados os 4 sinais anteriores e cada um");
disp("# deles será 'elevado ao quadrado' (amostra por");
disp("# amostra).");
disp("# ");
disp("Pressione qualquer tecla...");
input("");

disp("# ");
disp("# Relembrando o primeiro exemplo...");
plotnplay(SR,sinal1);
disp("# ... e elevando o sinal ao quadrado:");
plotnplay(SR,sinal1.*sinal1);
disp("# Relembrando o segundo exemplo...");
plotnplay(SR,sinal2);
disp("# ... e elevando o sinal ao quadrado:");
plotnplay(SR,sinal2.*sinal2);
disp("# Relembrando o terceiro exemplo...");
plotnplay(SR,sinal3);
disp("# ... e elevando o sinal ao quadrado:");
plotnplay(SR,sinal3.*sinal3);
disp("# Relembrando o quarto exemplo...");
plotnplay(SR,sinal4);
disp("# ... e elevando o sinal ao quadrado:");
plotnplay(SR,sinal4.*sinal4);



disp("# ");
disp("# Nada impede que multipliquemos quaisquer dois");
disp("# sinais, por exemplo, uma senóide em 1000Hz pela");
disp("# música do quarto exemplo:");
disp("#");
disp("Pressione qualquer tecla...");
input("");

# para criar uma senóide com a mesma duração do sinal 4, primeiro
# calculamos o número de amostras dele:
n=length(sinal4);
# e passamos como terceiro argumento o valor n/SR
# duração em segundos = no. de amostras do sinal / no. de amostras por segundo
plotnplay(SR,sinetone(1000,SR,n/SR,1).*sinal4);





disp("# ");
disp("# O próximo exemplo ilustra uma aplicação do produto");
disp("# não dos sinais em função do tempo, mas de seus");
disp("# espectros. São multiplicados os espectros de uma");
disp("# versão do exemplo 3 (senóides com frequências");
disp("# aleatórias) com a música do exemplo 4. O espectro");
disp("# do exemplo 3 é 'suavizado' tomando-se as médias em");
disp("# blocos de 200 pontos. Não se preocupem em entender");
disp("# a fundo este exemplo, apenas concentrem-se nas");
disp("# características do som resultante.");
disp("#");
disp("Pressione qualquer tecla...");
input("");

# s é recriado com as frequências aleatórias do exemplo 3,
# porém com a duração do exemplo 4
n=length(sinal4);
s = sinetone(freq(1),SR,n/SR,0.3)+sinetone(freq(2),SR,n/SR,0.3)+sinetone(freq(3),SR,n/SR,0.3);
disp("# ");
disp("# Relembrando o exemplo 3...");
disp("# ");
plotnplay(SR,sinal3);
# x é o espectro (complexo) de s, e xmed é a versão 'suavizada' do mesmo.
# observe que usamos a mesma função média que definimos antes...
x = fft(s);
xmed = media(x,200);
# y é o espectro do sinal4, e ymed é apenas uma versão recortada dele
# com o mesmo tamanho de xmed (n-200+1), para compatibilizar o produto
# que iremos fazer.
y = fft(sinal4);
ymed = y(1:length(xmed));
# z é obtido fazendo-se o produto ponto-a-ponto dos espectros xmed e ymed,
# aplicando-se a transformada inversa de Fourier, e convertendo os valores
# para números reais (eles são "quase" reais, mas possuem alguns valores
# minúsculos nas componentes imaginárias).
z = real(ifft(xmed.*ymed));
# z deve ser "normalizado" para a faixa [-1...+1]
z = z/max(abs(z));
disp("# ");
disp("# E a transformação do exemplo 4 através do exemplo 3");
disp("# ");
plotnplay(SR,z);








disp("# ");
disp("# O exemplo a seguir ilustra uma técnica de síntese");
disp("# conhecida como");
disp("# síntese por formantes. Nesta versão simplificada");
disp("# são calculadas somas de senóides como nos exemplos");
disp("# 2 e 3, porém com frequências cuidadosamente");
disp("# escolhidas. O resultado é uma mudança de coloração");
disp("# (timbre) em cada um dos sons, porém junto com o");
disp("# timbre podemos perceber alguma informação");
disp("# adicional.");
disp("# ");
disp("Pressione qualquer tecla...");
input("");

# a matriz freq é inicializada com algumas frequências
freq = [ 730 1090 2440;
530 1840 2480;
270 2290 3010;
570 840 2410;
300 870 2240];
# os sons terão todos a fundamental de 110Hz
fund = 110;
# para cada linha da matriz, será construido um som
sinal=[];
for i=1:5
  # constituido pela fundamental com amplitude 1/2
  s = sinetone(fund,SR,1,0.5);
  for j = 1:3
    # e por harmônicos próximos aos valores da tabela acima.
    # observe que esta fórmula corresponde a arredondar os
    # valores da tabela à série harmônica 110, 220, 330, etc.
    f = round(freq(i,j)/fund)*fund;
    # s é misturado ao harmônico correspondente com amplitude 1/6
    s += sinetone(f,SR,1,0.133);
  endfor
  sinal = [sinal; s];
endfor

plotnplay(SR,sinal);