Reply by Jeff Brower ●March 31, 20042004-03-31
Reshma-

I have posted your question on the MATLAB Yahoo Group.

-Jeff Reshma Jacob wrote:
>
> Hello Sir,
> i am trying to implement the yule walk filter in C++.i have the IEEE paper"the
> modified yulewalker method of ARMA spectral estimation" on which the yulewalk
> filter is based.But the standard of the paper is beyond my understanding.So i
am
> relying mainly on the matlab code of yulewalk function in the signal
processing
> tool box.I have some doubts on why some steps are being done....If u can clear
my
> doubts it would be of gr8 help..
> i am giving the code below....i have indicated where i ahve doubts in the code
with
> '???????'......Sir can u give me anytips on wat to read for a greater
understanding
> on this topic..i have gone thru the text book by proachis where i ahve
referred to
> the chapter on power spectral estimation for info on yulewalk....but i found
it
> insufficient to understand how the filter actually works....
> if (nargin < 3 | nargin > 5)
> error('Wrong number of input parameters.')
> end
> if (nargin > 3)
> if round(2 .^ round(log(npt)/log(2))) ~= npt
> % NPT is not an even power of two
> npt = round(2.^ceil(log(npt)/log(2)));
> end
> end
> if (nargin < 4)
> npt = 512;
> end
> if (nargin < 5)
> lap = fix(npt/25);
> end
>
> [mf,nf] = size(ff);
> [mm,nn] = size(aa);
> if mm ~= mf | nn ~= nf
> error('You must specify the same number of frequencies and amplitudes.')
> end
> nbrk = max(mf,nf);
> if mf < nf
> ff = ff';
> aa = aa';
> end
> if abs(ff(1)) > eps | abs(ff(nbrk) - 1) > eps
> error('The first frequency must be 0 and the last 1.')
> end
> % interpolate breakpoints onto large grid
> npt = npt + 1; % For [dc 1 2 ... nyquist].
> Ht = zeros(1,npt);
> nint=nbrk-1;
> df = diff(ff');
> if (any(df < 0))
> error('Frequencies must be non-decreasing.')
> end
> nb = 1;
> Ht(1)(1);
> for i=1:nint
> if df(i) == 0
> nb = nb - lap/2;
> ne = nb + lap;
> else
> ne = fix(ff(i+1)*npt);
> end
> if (nb < 0 | ne > npt)
> error('Too abrupt an amplitude change near end of frequency interval.')
> end
> j=nb:ne;
> if ne == nb
> inc = 0;
> else
> inc = (j-nb)/(ne-nb);
> end
> Ht(nb:ne) = inc*aa(i+1) + (1 - inc)*aa(i);??????????????????why is this
> done...is to give more weight to a(i) for the earlier points and greater
weight for
> a(i+1) in later points...if it is so..why is this???
> ; nb = ne + 1;
> end
> Ht = [Ht Ht(npt-1:-1:2)];?????????????????????????????wats the purpose of
making it
> to a 1024 point array....how does this affect the result
> n = length(Ht);
> n2 = fix((n+1)/2);
> nb = na;
> nr = 4*na;??????????????wat implication does multiplication of order by 4
have...in
> my program i have used na ie order "...
> nt = 0:1:nr-1;
> % compute correlation function of magnitude squared response
> R = real(ifft(Ht .* Ht));???????inverse fourier transform of power spectral
density
> gives correlation....is this wats goin on in this step...
> R = R(1:nr).*(0.54+0.46*cos(pi*nt/(nr-1))); ??????????? I:NR points are
> chosen...ie 4 times the order ...why is this.... % pick NR correlations
> % Form window to be used in extracting the right "wing" of two-sided
> % covariance sequence.
> Rwindow = [1/2 ones(1,n2-1) zeros(1,n-n2)];
> A = polystab(denf(R,na)); % compute denominator
> Qh = numf([R(1)/2,R(2:nr)],A,na); % compute additive decomposition
> Ss = 2*real(freqz(Qh,A,n,'whole'))'; % compute impulse response
> hh = ifft(exp(fft(Rwindow.*ifft(log(Ss)))));
> B = real(numf(hh(1:nr),A,nb));
> function b = numf(h,a,nb)
> %NUMF Find numerator B given impulse-response h of B/A and denominator A
> % NB is the numerator order. This function is used by YULEWALK.
>
> % Was Revision: 1.3, Date: 1994/01/25 17:59:33
> nh = max(size(h));
> impr = filter(1,a,[1 zeros(1,nh-1)]);
> b = h/toeplitz(impr,[1 zeros(1,nb)])';
> function A = denf(R,na)
> %DENF Compute denominator from covariances.
> % A = DENF(R,NA) computes order NA denominator A from covariances
> % R(0)...R(nr) using the Modified Yule-Walker method.
> % This function is used by YULEWALK.
> % Was Revision: 1.5, Date: 1994/01/25 17:59:00
> nr = max(size(R));
> Rm = toeplitz(R(na+1:nr-1),R(na+1:-1:2));??????????wtas the logic behind
choosing a
> toeplitz matrix of this order
> Rhs = - R(na+2:nr);
> A = [1 Rhs/Rm'];