Matlab - Transfer Function + modeling

Author

Transfer Function + modeling

Zaki N.

2005-04-06, 12:53 pm

Hi everybody,

I'm working on the modeling of an acoustic room TF (Transfer
Function) with an ARMA model, so I made a program which uses a
least-mean squares method based on this work:

www.kecl.ntt.co.jp/icl/signal/makin...4sap320-328.pdf

To simplify the simulation, I've put:
-In the estimated TF (i.e. 'hj'): a 10* random vector.
-In the measured TF (i.e. 'h'): hj + noise.

When I run my program; I always obtain the same results as the
"prony" command, but never like the "stmcb" command, why?!

Can you please optimize/correct my program?

Please try my program with different (P, Q, N) values.

-Below my program.

Thanks in advance!

Zaki N.

2005-04-06, 12:53 pm

% --------------------------------------------------------
% --------------------------------------------------------
% Ordinary least-mean squares method
% ARMA model identification

clc, clear all, close all

N = 5; % Impulse response''s order''s
P = 6; % Denominator order''s
Q = 5; % Numerator order''s

% P >= Q

% ------------------------------------------- Fj
% Fj size's: (N+ P+ 1, P)

Fj = sparse (zeros (N+ P+ 1, P));

hj = sparse (round (10* (rand ((N+ 1), 1))));
% hj: estimated transfer function

for i = 1:P
Fj (:, i) = [zeros(i, 1); hj; zeros(P- i, 1)];
end
% ------------------------------------------- end of Fj

% ------------------------------------------- h
% h size's: ((N+ P+ 1), 1)
% h: mesured transfer function

h = zeros (N+ P+ 1, 1);

h (1: N+ 1, 1) = hj (:, 1) + randn (1, 1);

h = sparse (h);
% ------------------------------------------- end of h

% ------------------------------------------- D
% D size's: (N+ P+ 1, Q+ 1)

D = speye (N+ P+ 1, Q+ 1);
% ------------------------------------------- end of D

% ------------------------------------------- F
% F size's: (N+ P+ 1, P+ Q+ 1)
% F contains the estimated transfer function: Fj

F = zeros (N+ P+ 1, P+ Q+ 1);

ki = 0; kj = 0;
while (ki < (N+ P+ 1))
for j = P+ 1 : (P+ 1+ Q)
for i = 1 : (N+ P+ 1)
F (i+ ki, j+ kj) = D (i, j- P);
end
end
kj = kj + (Q+ 1);
ki = ki + (N+ P+ 1);
end

ki = 0;
while (ki < (N+ P+ 1))
for j = 1:P
for i = 1 : (N+ P+ 1)
F (i+ ki, j) = Fj (i, j);
end
end
ki = ki + (N+ P+ 1);
end

F = sparse (F);
% ------------------------------------------- end of F

% ------------------------------------------- X
% X size's: (P+ Q+ 1, 1)
% X: the unknown (numerator + denominator)

X = sparse ((F.' * F)\ (F.' * h));

for i = 1: P
a (i, 1) = -X (i, 1);
disp (['a(', num2str(i),') = ', num2str(a (i, 1))])
end
a = sparse (a);
% a: the denominator

in = 1; kj = 0;
for i = 1: (Q+ 1)
b (in, 1) = X (in+ P, 1);
disp (['b(', num2str(kj),') = ', num2str(b (in, 1))])
kj = kj + 1;
in = in + 1;
end
b = sparse (b);
% b: the numerator
% ------------------------------------------- end of X

% ------------------------------------------- e
% e size's: (N+ P+ 1, 1)
% e: the error vector

e = h - (F * X);
e = sparse (e);

J = e.' * e;
% ------------------------------------------- end of e

% ------------------------------------------- ARMA model
hj = full (hj);
a_j = full (a.');
b_j = full (b.');

a_j_olms = [1, a_j]; % a(0) = 1
b_j_olms (1: Q+ 1) = b_j (1: Q+ 1);
figure; freqz (b_j_olms, a_j_olms, 256)
title ('Ordinary least-mean squares method')
% 'my' method

[b_j_mod1, a_j_mod1] = prony (hj, Q, P);
figure; freqz (b_j_mod1, a_j_mod1, 256)
title ('Prony''s method')
% prony command

[b_j_mod2, a_j_mod2] = stmcb (hj, Q, P, 5);
figure; freqz (b_j_mod2, a_j_mod2, 256)
title ('Steiglitz-McBride''s iteration')
% stmcb command
% ------------------------------------------- end of ARMA

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

Zaki N.

2005-04-06, 12:53 pm

UP UP UP !

Anders Björk

2005-04-06, 12:53 pm

Is it the same or similar? Look up examples for known models/data for use
with your ARMA model. Hayes book on Digital Signal Processing has examples
if I rembember right.

BR
Anders

"Zaki N." <zaki_nabil@caramail.com> skrev i meddelandet
news:ef00f1e.-1@webx.raydaftYaTP...
> Hi everybody,
>
>
> I'm working on the modeling of an acoustic room TF (Transfer
> Function) with an ARMA model, so I made a program which uses a
> least-mean squares method based on this work:
>
> www.kecl.ntt.co.jp/icl/signal/makin...4sap320-328.pdf
>
>
>
> To simplify the simulation, I've put:
> -In the estimated TF (i.e. 'hj'): a 10* random vector.
> -In the measured TF (i.e. 'h'): hj + noise.
>
>
> When I run my program; I always obtain the same results as the
> "prony" command, but never like the "stmcb" command, why?!
>
> Can you please optimize/correct my program?
>
> Please try my program with different (P, Q, N) values.
>
> -Below my program.
>
> Thanks in advance!

Zaki N.

2005-04-06, 12:54 pm

Thanks for your answer Anders!

-With 'my' program I obtain similar results as the "prony" command;
but the main question is:

Why (with my Matlab 6.5) "stmcb" and "prony" never provide similar
results?!

I'm awaiting same results, when I use them in this form:
prony (hj, Q, P)
stmcb (hj, Q, P, 5)

:-(

Anders Björk

2005-04-06, 12:54 pm

Well look these methods up in good book on digital signal processing! I'm
sure a library near you have some books.
BR
Anders

"Zaki N." <zaki_nabil@caramail.com> skrev i meddelandet
news:ef00f1e.3@webx.raydaftYaTP...
> Thanks for your answer Anders!
>
> -With 'my' program I obtain similar results as the "prony" command;
> but the main question is:
>
> Why (with my Matlab 6.5) "stmcb" and "prony" never provide similar
> results?!
>
> I'm awaiting same results, when I use them in this form:
> prony (hj, Q, P)
> stmcb (hj, Q, P, 5)
>
> :-(

Sponsored Links