Sign in

Search compdsp

comp.dsp by Keywords

Sponsor

Industry's highest performing at the lowest power DSPs now as low as $5.00*
Start development today!
*volume pricing for 10ku

Discussion Groups

Free Online Books

Discussion Groups | Comp.DSP | Rootmusic

There are 2 messages in this thread.

You are currently looking at messages 0 to 2.

Rootmusic - Derrick Rea - 2004-06-04 11:59:00

This posting is a continuation of an earlier thread I was involved in.
I wanted to continue with the old thread but found that re-posting
wasn't possible. I guess the newsgroup world moves faster than me. I
hope this post makes sense outside the original. The original thread
can be found at -;

http://groups.google.co.uk/groups?hl=en&lr=&ie=UTF-8&threadm=5d5e74f5.0312100940.2ab4f7f6%40posting.google.c
om&rnum=1&prev=/groups%3Fq%3Drootmusic%26hl%3Den%26lr%3D%26ie%3DUTF-8%26selm%3D5d5e74f5.0312100940.2ab4f7f6%2540
posting.google.com%26rnum%3D1

I am working on analysis blood velocity signals using the rootmusic
algorithm. My work is Matlab based and I was having problems trying to
get the function to work properly. I have, over the last few months,
resolved most of the problems I asked for advice on through the other
posts in the thread above, and hope this may be of help to anyone else
working with pmusic(), rootmusic().
I have a solution to a problem that I think is rooted in Matlab code
and a better understanding of why the autocorrelation matrix size
affects the frequency resolution of the music pseudospectrum and
rootmusic() z-plane zero location.
In other rootmusic threads there have been suggestions that there is a
better way to solve the frequency resolution problem I talk about
before. If anyone knows how, can you please get back to me.
Rune, you mentioned this problem in a previous posting (below). If you
read this could you explain further. Thanks to you for your earlier
help.


I found MATLAB rootmusic(), pmusic() could correctly calculate the
eigenvectors and eigenvalues and was reasonable at estimating the
correct frequency locations of exponentials either through
identification of zeros in the z-plane or peaks in the pseudospectrum
but that a mistake in code was made in the rootmusic() subfunction
compute_power(). It stated the since the signal was being modelled
using real sinusoids, rootmusic() needed to use only half the
eigenvectors/eigenvalues to calculate the corresponding powers.
Matlab's doesn't then choose the correct eigenvectors/eigenvalues that
are representative of the complete real sinusoidal model. If the model
order was defined to be 40, rootmusic() chooses the 20 largest
eigenvalues/eigenvectors to resolve the powers through the matrix
equation below.

 p
sum Pk|Vi(e(jwk))| ^2 = signaleigenvalue inoise variance ;i=1,2..p
k=1

Equation from Monsoon Hayes: Statistical Digital Signal Processing and
Modelling. Equation 8.159 pg 460.

This is wrong. In symmetrical sinusoidal spectra each harmonic is made
up of a power at a +ve frequency and a power at its mirror frequency
on the ve axis. The two powers will be identical and through this
eigenspectrum analysis will have associated eigenvalues that are very
close in value. It is not possible to discriminate between negative
and positive halves of the spectrum simply by picking out strongest
eigenvalues.

I've taken out the line of code that performs this discrimination and
there are no more spurious high negative and positive powers.
There may be an application where this assumption is valid but not in
mine. The change to the MATLAB code required to produce decent spectra
without spurious powers is given below. If anyone thinks I'm talking
rubbish please get back to me.

>%-----------------------------------------------------------------------------------------------
>function powers = >compute_power(signal_eigenvects,eigenvals,w_i,p_eff,sigma_w,xIsReal)

>% For real sinusoids, the system of eqs. has half the number of
unknowns
>if xIsReal,
>   w_i = reshape(w_i,2,length(w_i)./2);
>   w_i = w_i(1,:); % Use only the positive freqs.
>   w_i = w_i(:);
>   p_eff = p_eff./2;  %  Remove this line to correct the error in this function as 
        % described above.
>end

>% Form the A matrix
>if length(w_i) == 1,
>   % FREQZ does not compute the gain at a single frequency, handle this separately
>   A = polyval(signal_eigenvects,exp(i*w_i));
>else
>   for n = 1:p_eff,
>      A(:,n) = freqz(signal_eigenvects(:,n),1,w_i);
>   end
>end

>A = abs(A.').^2;

>% Form the b vector
>b = eigenvals(1:p_eff) - sigma_w;
>w_i_degrees = w_i * (200./(2*pi)); %added for debugging.

>% The powers are simply the solution to the set of eqs.
>powers = lsqnonneg(A,b); %forces non negative powers.
>%powers = A\b;

The second part of my problem when trying to apply rootmusic() was
what size of autocorrelation function to use. This is my posting and
Rune's answer.

http://groups.google.co.uk/groups?hl=en&lr=&ie=UTF-8&selm=f56893ae.0312151753.7f864a11%40posting.google.com

>> Increasing the size of the
>> autocorrelation matrix should only increase the number of noise
>> eigenvector and so help locate the signal frequencies.
>> 
>> I investigated the relationship between the correlation matrix size
>> and MUSICs ability to locate my 12 complex frequencies and I found
>> that the likelyhood of finding the frequencies increases as the
size
>> of the autocorrelation matrix is increased. Is generally true or is
>> there a limit of lag when the estimation of the autocorrelation
>> sequence becomes so inaccurate that MUSIC finds complex
exponentials
>> that aren't there. Why did you advice before that I would need to
work
>> with an order of approx 20 for 12 complex exponentials buried in
>> noise. Isn't bigger better?

>No, it isn't for a number of reasons. 

>First, you don't want to compute eigen vector decompositions of
>larger
>matrixes than absolutely nevessary. The computational task of
>decompsing
>a 200 x 200 matrix is significantly higher than decomposing a 20 x
>20
>matrix. It's a question of both time and numerical accuracy.

>Second, you want to pay some attention to the individual entries in
>the
>autocovariance matrix. If you compute the aotocovariance sequence
>the
>"usual" way, you implement something like (the process is assumed to
>be
>zero mean)
>
>                  1      N-1-|k|
>   Rxx(k) =    -------    sum    x(n)x^*(n+|k|)
>                N-|k|     n=0
>
>which means that only N-|k| cross-terms are summed to produce the 
>autocovariance coefficient of lag |k|. All available data points are 
>used to compute the autocovariance coefficient of lag 0 while only 
>x(0) and x(N-1) are used to compute autocovariance coefficients of 
>lags and -N+1 and N-1. Which, in turn, means that there is different 
>bias to estimates of coefficients of different lags. The basic MUSIC
>recipe does not address the question of biased estimators of the 
>signal autocovariance, so it is reasonable to believe that MUSIC 
>assumes the estimator to be unbiased.

>It is possible to make unbiased estimators for the signal
>autocovariance
>matrix for orders up to N/2. This unbiasedness comes at the expence
>of
>higher variance of the individual terms. 

>In summary: 

>- You want a low order of the autocovariance matrix for
>computational
>  reasons
>- You want an unbiased autocovariance estimator
>- You want as little variance of the autocovariance estimator as
>  possible. 
>- You still need the basic requirements, that relate the number of
>  sinusoidals present to the order of the autocovariance matrix, to 
>  hold, with a little slack for the unexpected.

>Finding both an estimator and an order that meets all these criteria 
>isn't always easy, and I usually say that these things have more to 
>do with black art and voodoo than anything else. My experience is 
>that using an order P such that

>  1.5*M < P < 2*M

>where M is the true (maximum) number of sinusoidals present, usually 
>works quite well. There are, of course, exceptions that go both
>ways,
>but the above works well as a general rule of thumb.

I was asking for help on why the size of the autocorrelation matrix
affects the frequency resolution of the rootmusic() algorithm. I have
discovered that the reason that the size of the autocorrelation matrix
has been so important in my application is that MATLAB calculates the
pseudospectrum as using the sum of ffts of the noise eigenvectors.
This obviously incurs the windowing limitations apparent in the
periodogram.

I found that the autocorrelation matrix sizes indicated above ie
1.5*M < P < 2*M weren't large enough to resolve the harmonic
frequencies in my signal. Instead I needed to supply a correlation
matrix as long as one complete cycle of the fundamental frequency in
my signal. One period of my blood velocity signal was 220 samples long
and so I needed to evaluate an autocorrelation matrix of 220x220 to
resolve frequencies separated by 1 / 220 * sample interval.
This dependency of the pseudoperiodogram on noise eigenfilter fft
seems wrong.
Rune recognises the problem that using the DFT introduces this
limitation with frequency resolution and questions that this is a
proper method in the posting below. Is there another approach to
calculating the pseudospectrum from the music equation below. I
thought about generating a test matrix [1 e(-jwn t)  e(-j(M-1)wn t)]
for all values of wn in my pseudospectrum but this appears to be
exactly the same as taking the sum of fft of the noise eigenvector.
Can anyone help. Rune are you still out there.

              1
    P(f) = ---------------
            M
            sum |Vi^H*e(f)|
           i=p+1

Vi  eigenvectors, e(f)  complex exponential test vector, [1
e(-j2pift)  e(-j2(M-1)pift)]

http://groups.google.co.uk/groups?hl=en&lr=&ie=UTF-8&selm=f56893ae.0304260209.5b4abd31%40posting.google.com

>f...@posting.google.com

>k...@yahoo.com (Jay) wrote in message 

>> These are sorted
>> and the spectrum is given by
>> reciprocal of the fft of the Sum of the Eigen Vectors from the
noise
>> subspace.

>Eh... that was a way of formulating the MUSIC _pseudo_ spectrum that 
>was new to me. At first glance, this formulation appears to limit the
>resolution of MUSIC to the DFT resolution of a sequence of the same 
>length as the eigen vectors. Are you sure the DFT is a proper method 
>for computing the pseudo spectrum?

Derrick

______________________________
New DSP Code Snippets Section now Live. Learn more about the reward program for contributors here.

Re: Rootmusic - Rune Allnor - 2004-06-05 09:46:00

d...@hotmail.com (Derrick Rea) wrote in message
news:<5...@posting.google.com>...
> This posting is a continuation of an earlier thread I was involved in.
> I wanted to continue with the old thread but found that re-posting
> wasn't possible. I guess the newsgroup world moves faster than me. I
> hope this post makes sense outside the original. The original thread
> can be found at -;
> 
>
http://groups.google.co.uk/groups?hl=en&lr=&ie=UTF-8&threadm=5d5e74f5.0312100940.2ab4f7f6%40posting.google.c
om&rnum=1&prev=/groups%3Fq%3Drootmusic%26hl%3Den%26lr%3D%26ie%3DUTF-8%26selm%3D5d5e74f5.0312100940.2ab4f7f6%2540
posting.google.com%26rnum%3D1
> 
> I am working on analysis blood velocity signals using the rootmusic
> algorithm. My work is Matlab based and I was having problems trying to
> get the function to work properly. I have, over the last few months,
> resolved most of the problems I asked for advice on through the other
> posts in the thread above, and hope this may be of help to anyone else
> working with pmusic(), rootmusic().
> I have a solution to a problem that I think is rooted in Matlab code
> and a better understanding of why the autocorrelation matrix size
> affects the frequency resolution of the music pseudospectrum and
> rootmusic() z-plane zero location.
> In other rootmusic threads there have been suggestions that there is a
> better way to solve the frequency resolution problem I talk about
> before. If anyone knows how, can you please get back to me.
> Rune, you mentioned this problem in a previous posting (below). If you
> read this could you explain further. Thanks to you for your earlier
> help.

I based my work on the Tufts/Kumaresan version of Prony's method.
It is implemented much the same way as Root MUSIC, but builds a 
somewhat different polynomial that is more directly linked to the 
signal, and that also is easier from a pure numerical point of view. 
The TK polynomial is half the order of the Root MUSIC polynomial, 
which can be significant for higher orders.

Check out 

Tufts and Kumaresan: Estimation of frequencies of multiple sinusoids:
   Making Linear Prediction perform like Maximum Likelihood.
   Proc. IEEE, vol 70 no 9,September 1982, p 975.     

Another good reference would be 

Therrien: Discrete Random Signals and Statistical Signal Processing
         Prentice-Hall, 1992

if you can find a copy. Unfortunately, this book has been out of stock 
for almost a decade.

> I was asking for help on why the size of the autocorrelation matrix
> affects the frequency resolution of the rootmusic() algorithm. I have
> discovered that the reason that the size of the autocorrelation matrix
> has been so important in my application is that MATLAB calculates the
> pseudospectrum as using the sum of ffts of the noise eigenvectors.
> This obviously incurs the windowing limitations apparent in the
> periodogram.
> 
> I found that the autocorrelation matrix sizes indicated above ie
> 1.5*M < P < 2*M weren't large enough to resolve the harmonic
> frequencies in my signal. Instead I needed to supply a correlation
> matrix as long as one complete cycle of the fundamental frequency in
> my signal. One period of my blood velocity signal was 220 samples long
> and so I needed to evaluate an autocorrelation matrix of 220x220 to
> resolve frequencies separated by 1 / 220 * sample interval.
> This dependency of the pseudoperiodogram on noise eigenfilter fft
> seems wrong.
> Rune recognises the problem that using the DFT introduces this
> limitation with frequency resolution and questions that this is a
> proper method in the posting below. Is there another approach to
> calculating the pseudospectrum from the music equation below. I
> thought about generating a test matrix [1 e(-jwn t) ? e(-j(M-1)wn t)]
> for all values of wn in my pseudospectrum but this appears to be
> exactly the same as taking the sum of fft of the noise eigenvector.
> Can anyone help. Rune are you still out there.

Yep I am. I tried to make some precautinary remarks that I might be 
wrong, as can be seen in the quotation of my post below.It turned out 
that I was indeed wrong. The DFT is perfectly usable for finding close 
*zeros* in a spectrum. The resolution limitation I was thinking of, 
applies to finding spectrum *lines*, which is a completely different 
cup of tea. Mea Culpa.

I've made a short matlab script that demonstrates how to find 
close zeros from short polynomials:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
R=[1,1];                 % Magnitude of zero
TH=[0.14,0.2]*2*pi;      % Phase of zero

[re,im]=pol2cart(TH,R);
p=re+i*im;

P=poly(p);               % P is complex-valued polynomial of order 2

r1=roots(P);
[th,r]=cart2pol(real(r1),imag(r1));

% Consistency check: See if P really has its roots where I specified 
% them to be

clf
subplot(211)
polar(TH,R,'or')         % Plot designed roots
hold on
polar(th,r,'xb')         % Plot roots of P

% Red circles and blue x'es overlap: Design OK.

% Zero-pad P to total length of 1024 and compute then plot sepctrum.

N=1024;
fv=[0:N-1]/N;
S=log(abs(fft(P,N)));

subplot(212)
plot(fv,S)
% Spectrum has zeros at 0.14*2pi and 0.2*2pi, as specified above
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
>               1
>     P(f) = ---------------
>             M
>             sum |Vi^H*e(f)|
>            i=p+1
> 
> Vi ? eigenvectors, e(f) ? complex exponential test vector, [1
> e(-j2pift) ? e(-j2(M-1)pift)]
> 
>
http://groups.google.co.uk/groups?hl=en&lr=&ie=UTF-8&selm=f56893ae.0304260209.5b4abd31%40posting.google.com
> 
> >f...@posting.google.com
>  
> >k...@yahoo.com (Jay) wrote in message 
>  
> >> These are sorted
> >> and the spectrum is given by
> >> reciprocal of the fft of the Sum of the Eigen Vectors from the
>  noise
> >> subspace.
>  
> >Eh... that was a way of formulating the MUSIC _pseudo_ spectrum that 
> >was new to me. At first glance, this formulation appears to limit the
> >resolution of MUSIC to the DFT resolution of a sequence of the same 
> >length as the eigen vectors. Are you sure the DFT is a proper method 
> >for computing the pseudo spectrum?
> 
> Derrick

Rune