## Contents

% Blog for linear prediction lecture

format shortg


## input the data from a .wav file and select some data to analyze

[Y1,FS]=audioread('BeeGees.wav');
NBITS = 16;
% transform the signal from (-1,1) into (-2^(NBITS-1), 2^(NBITS-1))
Y = Y1*2^(NBITS-1);
%select some data to analyze
% select a small segment containing only four seconds of audio
Y = Y1*2^(NBITS-1);
Tbase = 1000000+2*44100;
figure(1),clf
plot(Y(Tbase+(1:4*44100)))
Ya = Y(Tbase+(1:4*44100));
% play in audioplayer only the selected segment
player = audioplayer(Ya/2^NBITS,FS);
playblocking(player);

% take segments of 80 ms each
segs = 1:(44100*0.08);
Lsegs = length(segs);
Ys = Ya(segs);
figure(2),plot(Ys)
title('Audio segment to analyze using L-D algorithm')
xlabel('time')
% Extract the mean
Ysm = Ys- mean(Ys);
% Estimate the autocorrelation function for the signal in the current block
for tau = 0:300
r(tau+1) = sum(Ysm(1:(Lsegs-tau)).* Ysm(tau+(1:(Lsegs-tau))))/Lsegs;
end
r = r(:); % be sure that r is a column vector
figure(3),plot(r)
grid
title('estimated autocorrelation function ')
xlabel('lag \tau')
ylabel('r(\tau)')


## Construct the optimal forward predictor for order m = 3

m = 3
% the autocorrelation matrix of order m
R = toeplitz(r(1:(m)))
% the vector r
r_vec = r(2:(m+1));
% the optimum prediction filter w
w = R\r_vec
% the variance of prediction error
P = r(1) - r_vec'*w
% the prediction error filter
a = [1; -w]
T1 = toeplitz(r(1:(m+1)))
% the augmented Wiener-Hopf equation is T1*a = [P; zeros(m,1)]
T1*a
% check the augmented Wiener-Hopf equation
T1*a - [P; zeros(m,1)]
% the augmented Wiener-Hopf equation for the backward predictor
% T1*a(end:-1:1) = [zeros(m,1); P]
% check the identity
T1*a(end:-1:1) - [zeros(m,1); P]

m =

3

R =

4.3064e+06    4.301e+06   4.2898e+06
4.301e+06   4.3064e+06    4.301e+06
4.2898e+06    4.301e+06   4.3064e+06

w =

1.507
-0.46253
-0.046476

P =

7710.4

a =

1
-1.507
0.46253
0.046476

T1 =

4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06
4.301e+06   4.3064e+06    4.301e+06   4.2898e+06
4.2898e+06    4.301e+06   4.3064e+06    4.301e+06
4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06

ans =

7710.4
1.4359e-09
-9.7727e-11
8.2114e-10

ans =

1.0186e-10
1.4359e-09
-9.7727e-11
8.2114e-10

ans =

9.3132e-10
0
1.8626e-09
0



## Check the Levinson-Durbin recursion from m= 3 to m = 4

all relevant variable for m= 3 have been computed above use subscript _n for the variables at m= 4 and compute directly the predictors from their definitions

r_vec_n = r(2:(m+2));
R_n = toeplitz(r(1:(m+1)))
w_n = R_n\r_vec_n
a_n = [1; -w_n]
% now use Levinson-Durbin recursion
% the definition of the crosscorrelation Delta
Delta = a'*r_vec_n(end:-1:1)
% the definition of the reflection coefficient
Gamma = -Delta/P
% the predictor a_n  compute by LD
psi = [a;0]+Gamma*[0; a(end:-1:1)]
% check that psi is equal to a_n
psi-a_n

R_n =

4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06
4.301e+06   4.3064e+06    4.301e+06   4.2898e+06
4.2898e+06    4.301e+06   4.3064e+06    4.301e+06
4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06

w_n =

1.4951
-0.58067
0.33841
-0.2554

a_n =

1
-1.4951
0.58067
-0.33841
0.2554

Delta =

-1969.2

Gamma =

0.2554

psi =

1
-1.4951
0.58067
-0.33841
0.2554

ans =

0
9.5257e-14
-3.5483e-13
4.8839e-13
-2.2893e-13



## The full Levinson-Durbin algorithm

Given data: autocorrelation sequence r

r = r(:);
% Initialize and run for m = 1:5
a = 1;
Delta = r(2);
P = r(1);
for m = 1:8
Gamma(m) = - Delta/P;
a = [a;0] + Gamma(m)*[0;a(end:-1:1)];
Delta = a'*r((m+2):-1:2);
P = P*(1-Gamma(m)^2);
end
% the final predictor is a
a'

ans =

Columns 1 through 6

1      -1.4146      0.51076     -0.22979    -0.022839     0.084494

Columns 7 through 9

0.057511    -0.060969     0.079223



## Compare the L-D solution to the direct one

The solution for m = 5 is a_1 as follows

R= toeplitz(r(1:(m)));
w = R\r(2:(m+1));
a_1 = [1; -w];
% check that a given by L-D and a_1 are the same
a_1-a

ans =

0
2.6157e-13
-4.5874e-13
1.9204e-13
6.347e-14
-2.6244e-13
2.5976e-13
5.9446e-14
-1.1474e-13



## Use the matlab built-in function levinson

[aM,PM,GammaM] = levinson(r,m);

% Check that Gamma and GammaM are the same
GammaM(:)-Gamma(:)
a_1-a

ans =

0
0
0
2.6423e-14
-6.1617e-14
-1.3617e-13
8.3183e-14
7.0249e-14

ans =

0
2.6157e-13
-4.5874e-13
1.9204e-13
6.347e-14
-2.6244e-13
2.5976e-13
5.9446e-14
-1.1474e-13



## Construct the optimal forward predictor for orders m \in {1,2,...,8}

Arrange all optimal predictors as rows below the main diagonal of a matrix U

U = zeros(8);
for m = 1:8
% the autocorrelation matrix of order m
R = toeplitz(r(1:(m)))
% the vector r
r_vec = r(2:(m+1));
% the optimum prediction filter w
w = R\r_vec
% the variance of prediction error
P = r(1) - r_vec'*w
a = [1; -w];
U(m+1,m+1+(0:-1:-m)) = a';
end
R = toeplitz(r(1:(m+1)))
% The resulting matrix U
U
% Check that U is diagonalizing the matrix R
U*R*U'

R =

4.3064e+06

w =

0.99874

P =

10805

R =

4.3064e+06    4.301e+06
4.301e+06   4.3064e+06

w =

1.5318
-0.53373

P =

7727.1

R =

4.3064e+06    4.301e+06   4.2898e+06
4.301e+06   4.3064e+06    4.301e+06
4.2898e+06    4.301e+06   4.3064e+06

w =

1.507
-0.46253
-0.046476

P =

7710.4

R =

4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06
4.301e+06   4.3064e+06    4.301e+06   4.2898e+06
4.2898e+06    4.301e+06   4.3064e+06    4.301e+06
4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06

w =

1.4951
-0.58067
0.33841
-0.2554

P =

7207.5

R =

4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06   4.2567e+06
4.301e+06   4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06
4.2898e+06    4.301e+06   4.3064e+06    4.301e+06   4.2898e+06
4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06    4.301e+06
4.2567e+06   4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06

w =

1.4421
-0.5104
0.21783
0.055062
-0.20765

P =

6896.7

R =

4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06   4.2567e+06   4.2336e+06
4.301e+06   4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06   4.2567e+06
4.2898e+06    4.301e+06   4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06
4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06    4.301e+06   4.2898e+06
4.2567e+06   4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06    4.301e+06
4.2336e+06   4.2567e+06   4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06

w =

1.4233
-0.50542
0.23751
0.0089484
-0.077357
-0.090348

P =

6840.4

R =

Columns 1 through 6

4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06   4.2567e+06   4.2336e+06
4.301e+06   4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06   4.2567e+06
4.2898e+06    4.301e+06   4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06
4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06    4.301e+06   4.2898e+06
4.2567e+06   4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06    4.301e+06
4.2336e+06   4.2567e+06   4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06
4.2064e+06   4.2336e+06   4.2567e+06   4.2752e+06   4.2898e+06    4.301e+06

Column 7

4.2064e+06
4.2336e+06
4.2567e+06
4.2752e+06
4.2898e+06
4.301e+06
4.3064e+06

w =

1.4187
-0.5094
0.23797
0.021162
-0.10335
-0.017155
-0.051424

P =

6822.3

R =

Columns 1 through 6

4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06   4.2567e+06   4.2336e+06
4.301e+06   4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06   4.2567e+06
4.2898e+06    4.301e+06   4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06
4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06    4.301e+06   4.2898e+06
4.2567e+06   4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06    4.301e+06
4.2336e+06   4.2567e+06   4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06
4.2064e+06   4.2336e+06   4.2567e+06   4.2752e+06   4.2898e+06    4.301e+06
4.1759e+06   4.2064e+06   4.2336e+06   4.2567e+06   4.2752e+06   4.2898e+06

Columns 7 through 8

4.2064e+06   4.1759e+06
4.2336e+06   4.2064e+06
4.2567e+06   4.2336e+06
4.2752e+06   4.2567e+06
4.2898e+06   4.2752e+06
4.301e+06   4.2898e+06
4.3064e+06    4.301e+06
4.301e+06   4.3064e+06

w =

1.4146
-0.51076
0.22979
0.022839
-0.084494
-0.057511
0.060969
-0.079223

P =

6779.5

R =

Columns 1 through 6

4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06   4.2567e+06   4.2336e+06
4.301e+06   4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06   4.2567e+06
4.2898e+06    4.301e+06   4.3064e+06    4.301e+06   4.2898e+06   4.2752e+06
4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06    4.301e+06   4.2898e+06
4.2567e+06   4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06    4.301e+06
4.2336e+06   4.2567e+06   4.2752e+06   4.2898e+06    4.301e+06   4.3064e+06
4.2064e+06   4.2336e+06   4.2567e+06   4.2752e+06   4.2898e+06    4.301e+06
4.1759e+06   4.2064e+06   4.2336e+06   4.2567e+06   4.2752e+06   4.2898e+06
4.1419e+06   4.1759e+06   4.2064e+06   4.2336e+06   4.2567e+06   4.2752e+06

Columns 7 through 9

4.2064e+06   4.1759e+06   4.1419e+06
4.2336e+06   4.2064e+06   4.1759e+06
4.2567e+06   4.2336e+06   4.2064e+06
4.2752e+06   4.2567e+06   4.2336e+06
4.2898e+06   4.2752e+06   4.2567e+06
4.301e+06   4.2898e+06   4.2752e+06
4.3064e+06    4.301e+06   4.2898e+06
4.301e+06   4.3064e+06    4.301e+06
4.2898e+06    4.301e+06   4.3064e+06

U =

Columns 1 through 6

0            0            0            0            0            0
-0.99874            1            0            0            0            0
0.53373      -1.5318            1            0            0            0
0.046476      0.46253       -1.507            1            0            0
0.2554     -0.33841      0.58067      -1.4951            1            0
0.20765    -0.055062     -0.21783       0.5104      -1.4421            1
0.090348     0.077357   -0.0089484     -0.23751      0.50542      -1.4233
0.051424     0.017155      0.10335    -0.021162     -0.23797       0.5094
0.079223    -0.060969     0.057511     0.084494    -0.022839     -0.22979

Columns 7 through 9

0            0            0
0            0            0
0            0            0
0            0            0
0            0            0
0            0            0
1            0            0
-1.4187            1            0
0.51076      -1.4146            1

ans =

Columns 1 through 6

0            0            0            0            0            0
0        10805     1.99e-09  -1.2296e-09   4.7294e-10   4.3656e-11
0   1.6549e-09       7727.1   2.2192e-09  -5.9663e-10   2.8194e-10
0  -9.1784e-10   2.0238e-09       7710.4   1.6753e-09  -9.3314e-10
0  -9.3015e-10   4.9707e-10   1.9059e-09       7207.5   3.0886e-09
0            0   4.6566e-10  -9.3458e-10   2.0155e-09       6896.7
0   -9.322e-10   7.2676e-10   3.8371e-11  -5.6069e-10   2.4751e-09
0   4.6522e-10   -9.468e-10   4.4914e-10  -3.1919e-12  -6.6895e-10
0   9.3132e-10  -4.9528e-10  -9.7273e-10   2.2562e-10   6.7717e-10

Columns 7 through 9

0            0            0
-8.5129e-10   5.8935e-10   1.1714e-09
6.4028e-10  -1.2187e-09  -4.3292e-10
3.0013e-10   2.0009e-10  -5.7844e-10
2.0009e-10  -9.0586e-10   3.2378e-10
2.261e-09  -1.1441e-09   1.5116e-09
6840.4    2.281e-09  -4.5111e-10
1.9504e-09       6822.3   2.3701e-09
-1.2619e-09   2.4493e-09       6779.5



## Check the dependency of the Error power, P, on the order of the system. Analyze the signal using "lpc" function

p = 1:300;
for ip = 1:length(p)
pc = p(ip);
[a,g(ip)] = lpc(Ysm,pc);
end
figure(4),plot(p,g,'or')