Re: loops vs vectorization

[Paul Kienzle]

On Jul 6, 2006, at 6:16 PM, Joan Picanyol i Puig wrote:

> 'k, so let's turn this into "I wish I'd know how to vectorize this".
>
> * Robert A.Macy <address@hidden> [20060706 01:11]:
> > for example, the simple act of replacing every item that is
> > less than negative with the value zero.
>
> Indexing with find should fit your needs.

This is easy:  x(x<0)=0

> I just found this in some code of mine:
>
>     %XXX this shouldn't be hard to vectorize
>     for lp = 1:lnechoes
>         ltmp(lp,llags(lp):ln_t + llags(lp) - 1) = ltmps(lp,:) * 
> lvalues(lp),
>     end

The right hand side is easy enough:
   diag(lvalues)*ltmps
The left hand side is harder:
   ltmp(sub2ind(size(ltmp),i,j))
where
   i=repmat([1:lnechoes],ln_t,1)
   j=repmat([1:ln_t]',1,lnechoes)+repmat(llags(:)',ln_t,1)


If what you want is really an attenuated echo, though, then you should 
be able to use filter.  E.g., the following is a lag 5 echo attenuated 
by 40%:

    x=zeros(1,40); x(1:2)=1; plot(filter(1,[1,0,0,0,-0.6],x))


> Any tips for this one? I have a harder one too, and this one is 
> killing me:
>
>     for i = 1:ln_t
>         [u,s,v] = svd(Rx(:,:,i)),
>         w = 1./diag(s + 1000*eps),
>         m = (abs(u.'*conj(a_s))).^2,
>         pow_density = 1 ./ (sum(diag(w)*m)),
>         k = find(pow_density == max(pow_density)),
>
>         %check for aliasing
>         if (length(k) ~= 1)
>             msg = sprintf('capon_doa_est found %d out of %d possible 
> DoAs\n',...
>                 k, length(theta_s))
>             warning(msg)
>             k = k(1)
>         end
>
>         ldoa(i) = theta_s(k),
>     end
>
> For those in signal processing, this is supposed to be a Capon DoA
> estimator.  The only reason I haven't rewritten it in C++ (besides the
> proverbial lack of time) is that I believe most time is spent on the
> svd() call, and thus there's not much to be gained. It is also svd()
> that I wouldn't know how to vectorize.

Check if repeating the svd doubles the amount of running time.  This 
will confirm that it is the bottleneck, and that the only improvement 
you will get is by changing your algorithm or using parallelism.

- Paul