% AR_PAROLE2.M
%
% P. Flandrin, 30 sept. 2003 - 4 mars 2005
%
% exemple de modŽlisation AR de signaux de parole
 
% ----- signal

% voise = wavread('Voise_nonvoise'); xx = voise'; Fs = 8192;
% phrase = wavread('Phrase'); xx = phrase'; Fs = 8192; 
% voyel_a = auread('voya'); xx = voyel_a'; Fs = 8192; 
% anna = wavread('annak_soleil'); xx = anna(300001:400000)'; Fs = 22300;

zoe = auread('zoe'); xx = zoe'; Fs = 44600; % load xZp.mat ; sound(xZp,Fs) % rejoue une synthese

% zoe = auread('zoe_8'); xx = zoe'; Fs = 16000; 
% aeiou = auread('aeiou'); xx = aeiou'; Fs = 16000;

N = length(xx);

figure(1)
clf
subplot(211)
plot(xx)
xlabel('temps')
set(gca,'XTick',[])
set(gca,'YTick',[])
axis([1 N 1.1*min(xx) 1.1*max(xx)])
hold on

pause

sound(xx,Fs)

[X1,Y] = ginput(1);
plot([fix(X1) fix(X1)],[1.1*min(xx) 1.1*max(xx)],'r')
[X2,Y] = ginput(1);
plot([fix(X2) fix(X2)],[1.1*min(xx) 1.1*max(xx)],'r')
x = xx(fix(X1):fix(X2));
plot(fix(X1):fix(X2),xx(fix(X1):fix(X2)),'r')
N = length(x);

hold off

subplot(212)

plot(x,'r')
axis([1 N 1.1*min(x) 1.1*max(x)])

pause

sound(x,Fs)

pause

% ----- identification

pmax = 20; % ordre max du mod�le estimŽ

[AA,P,FPE,popt] = AR_levinson(x,pmax);

figure(2)

subplot(311)
plot(log(P),'b+-')
hold on
plot(log(FPE),'r+-')
plot(popt,log(FPE(popt)),'r','LineWidth',6)
xlabel('ordre du modele')
ylabel('log(erreur)')
title('bleu = puissance; rouge = FPE')
hold off

subplot(312)
Nfft = 2048;
xf = fft(x,Nfft);
plot(linspace(0,.5*Fs,Nfft/2),log(abs(xf(1:Nfft/2))));
V = axis ;
axis([1 .5*Fs V(3) V(4)]) ;
ylabel('log(DSP_FFT)')
xlabel('frequence')

subplot(313)
Nf = 1000; % nombre de points du spectre
f = linspace(0,.5,Nf);
faff = linspace(0,.5*Fs,Nf);
z = exp(i*2*pi*f);

M = [];
M(1,:) = ones(1,Nf);
for mm = 2:pmax+1
	
	M(mm,:) = 1./z.^(mm-1);
	
end

Aest = [1 AA(popt,1:popt)];
Sest = 1./(abs(Aest*M(1:popt+1,:))).^2;

plot(faff,log(Sest),'r','LineWidth',2)
Ax = axis;
axis([0 .5*Fs Ax(3) Ax(4)])
xlabel('frequence')
ylabel('log(DSP_AR)')
title('pointilles = modele; trait plein = estimee') 

subplot(311)
hold on
[X,Y] = ginput(1);
pp = round(X(1));
% pp = 12;
plot(pp,log(P(pp)),'b','LineWidth',6)
hold off

subplot(313)
hold on
App = [1 AA(pp,1:pp)];
Spp = 1./(abs(App*M(1:pp+1,:))).^2;
plot(faff,log(Spp),'b','LineWidth',2)
hold off


%--------------

pe = 20 ;

Ae = [1 AA(pe,1:pe)];
ep = filter(Ae,1,x) ;

figure(3)

subplot(211)
plot(ep)
axis([1 N 1.1*min(ep) 1.1*max(ep)])
xlabel('temps')
title('erreur de prediction')

r = abs(xcorr(ep,'coeff')) ;
rr = r(length(ep):length(ep)+250-1) ;
subplot(212)
plot(rr,'+-')
xlabel('retard')
title('autocorrelation de l erreur de prediction')
[xdec,ydec] = ginput(2);
xdec1 = fix(xdec(1));
xdec2 = fix(xdec(2));
ind = find(abs(rr(xdec1:xdec2)) == max(abs(rr(xdec1:xdec2)))) + xdec1 - 1

xe = zeros(size(x)) ;
jnd = 1:ind:length(x) ;
xe(jnd) = 1 ;

xs = filter(1,[1 AA(pe,:)],xe) ;
xs = xs/max(abs(xs)) ;
xn = x/max(abs(x)) ;

figure(4)
subplot(211)
plot(xs)
axis([1 N 1.1*min(xs) 1.1*max(xs)])
xlabel('temps')
title('synthese')
subplot(212)
plot(xn)
axis([1 N 1.1*min(xn) 1.1*max(xn)])
xlabel('temps')
title('original')

pause
sound(xs,Fs)
pause
sound(x/max(abs(x)),Fs)

