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')