# Re: plz help

Started by March 31, 2004
 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'];