Mean Removal Experiments

June 29, 2007

This post is a follow up of the text posted yesterday (I’m afraid this will often be the case ;) ). It is related to radar, sonar and multi-channel signal processing, singular value decomposition, matched filtering, and more specifically to something that’s called clutter rejection.

Intro

There is another problem related to the one discussed yesterday. Consider again the function from example 3:

\sum_i A_i\exp{\left[-\jmath (x k_{x,i} + y k_y)\right]} = f(x)\cdot h(y),

where

f(x) = \sum_i A_i\exp{(-\jmath x k_{x,i})}

h(y) = \exp{(-\jmath y k_y)}

 

and

 

k_{x,i} \in [-K_x,K_x], \qquad K_x=\mbox{const.}

 

 

This gives a signal that can be decomposed into two singular vectors f(\mathbf{x}) and h(\mathbf{y}). Thus, one can write:

\mathbf{S} = \mathbf{f}\mathbf{h}

Suppose we know \mathbf{S} and k_y and we are trying to estimate A_i and k_{x,i}. Equation above implies that

 

 

\mathbf{\hat{f}} = \mathbf{S}\mathbf{h}^\dagger

 

 

\mathbf{\hat{f}} is an estimate of \mathbf{f} and ^\dagger denotes matrix pseudo inverse. This solution is called Wiener filter. For vectors, we have

 

 

\mathbf{f^\dagger} = \mathbf{h}'/n,

 

 

where ' denotes a conjugated transpose. n is the length of vector \mathbf{h}. In the signal processing literature, \mathbf{h}'/n is called matched filter. Further, suppose that k_y=0. In that case, \mathbf{h}=[1\, 1\, 1\, \dots]. Clearly, \mathbf{\hat{f}} can then be calculated as:

 

 

\mathbf{\hat{f}} = \frac{1}{n}\sum_j \{S\}_{i,j}.

 

 

This means simply a calculation of an average over columns of matrix \mathbf{S}.

 

 

Signal Separation

 

 

Now suppose that we have some other signals \mathbf{S_m} at unknown frequencies (k_{x,m},k_{y,m}) but k_{y,m}\neq k_y. These signals will be added to the signal from our previous section:

 

 

\mathbf{Z}=\mathbf{S}+\sum_m\mathbf{S_m}

 

 

One can then estimate signals \mathbf{S_m} like this:

 

 

\widehat{\sum_m\mathbf{S_m}}=\mathbf{Z}-\frac{1}{n}\mathbf{Z}\cdot\mathbf{h'}\mathbf{h}

 

 

Now, this leads to a least squares problem and that’s where the SVD comes into play; see the post from yesterday. I still need to figure this out…

 

 

Channel Mismatch

 

 

One can ask a question: “Why should I care about least squares if I know \mathbf{h}?” The answer is: You should. The reason is that the solution demonstrated above works only for the cases when we know \mathbf{h} (and therefore h(y)) precisely. This, however, isn’t the case in practice. In practice, function h(y) will be modulated by another function e(y), which is generally unknown. In radar, sonar and other multi-sensor applications, this function often models so-called channel mismatch. How does SVD help in this case? Well, what we know about the result r(y)=h(y)\cdot e(y) is that it is still a function of y only and thus it is orthogonal to function f(x). So, under the conditions discussed yesterday, we can still estimate functions f(x) and r(y) using SVD without really knowing them…

 

 

Da Code

 

 

Again, some code to play with. Here it is…

 

 


clear


cj=sqrt(-1);m=16;n=16;
x=linspace(-1,1,m)'*%pi;
y=linspace(-1,1,n)'*%pi;


s=zeros(m,n);
f=0;
Kx=7;ky=.1;
kx=linspace(-Kx,Kx,100);
h=exp(-cj*y'*ky);
for i=1:100,
g=exp(-cj*x*kx(i));
f=f+10*rand(1,1)*g;
end
s=f*h;


s1=zeros(m,n);
kx=-2;ky=2;
w=.1*exp(-cj*x*kx);
z=exp(-cj*y'*ky);
s1=s1+w*z;


kx=3;ky=-5;
w=.5*exp(-cj*x*kx);
z=exp(-cj*y'*ky);
s1=s1+w*z;


s=s+s1;


//f_=sum(s,2)/n;
f_=s*h'/n;
s=s-f_*h; //try commenting this out
//to see all signals


pad=96;
S=fftshift(fft2(fftshift([s,...
zeros(m,pad-n);zeros(pad-m,pad)])));


fig1=scf(1);
map=jetcolormap(256);
xset("colormap",map)
Matplot(log10(abs(S)+1e-12)*50)
Posted by doctuh4iah
Filed in Scientific Stuff
Leave a Comment »

Leave a Reply