Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.atnf.csiro.au/projects/askap/techdocs/eigen_filtering.pdf Дата изменения: Thu Oct 19 08:19:20 2000 Дата индексирования: Tue Oct 2 17:53:49 2012 Кодировка: Поисковые слова: m 101

Some Thoughts on Eigen Filtering
As a method of optimal filtering
Jon Bell CSIRO ATNF 1 Introduction The correlation matrix R of discrete time sampled data is the matrix formed when a vector of such data x(n) and its transpose are multiplied together:

И Й Й R=Й Й Й Й К

x0 x 0 x1 x 0 . . x N x0

x0 x1 x2 x2 . . .

. . x3 x . .

3

. x0 . x1 . . . xN

xN xN . . xN

Ы З З З З З З Ш

This matrix has eigen values i with corresponding eigen vectors qi such that (R - iI) qi = 0. Consider the following system
Signal u(n) Input x(n) Output y(n)


Noise v(n)

FIR Filter wn

Figure 1: Linear filtering to remove additive white Gaussian noise

The average power of the of the signal component of the output y(n) is P0 = average power of the noise component is N0 = 2wHw. We want to maximise filter) the output signal to noise ratio (SNR)0 = P0/N0 with the constraint th coefficients w0 of the impulse response of the optimal FIR filter in this case be: w0=q
max

wHRw, while the (find the optimal at wHw = 1. The can be shown to

where qmax is the eigen vector corresponding to the largest eigen value matrix R. These notes are based on "Adaptive Filter Theory" by Simon Haykin.

max

of the correlation

1


2 How well do such eigen filters work ? Figure 1 shows some synthetic test data u(n) to which random noise v(n) is added. The coefficients of the optimal filter derived as described above are shown in Figure 2, (with the additional step that the correlation matrix was averaged over 10 blocks of noisy data). Curiously this appears to be no more than a scaled version of the average of the signal over the 10 blocks of noisy data. Not surprising the spectral response of this filter looks like a noisy version of the spectral shape of the raw data as shown in Figure 3. Also shown in Figure 3 is spectral response of a smoothed version of this filter and digital filter formed using the Remez algorithm. These spectral responses suggest that the eigen filter is a rather bad filter when compared to the Remez filter. Figures 4 and 5 show the results of apply these two filters to the individual blocks of data, which were subsequently averaged. The result of the eigen filter is a very poor representation of the data. The result of the remez filter is much better (including a substantially lower RMS) when compared to a simple average of the data blocks, however, it does suffer from an increased DC content.

1 0.5 0 0 2 0 -2 0 1 0 -1 0

20

40

60

80

100

120

20

40

60

80

100

120

20

40

60

80

100

120

Figure 1: Test data. Top: raw data u(n). Middle: 1 block of noisy input data x(n). Bottom: Average over 10 blocks of noisy input data.

2


0 .3 0 .2 5 0 .2 0 .1 5 0 .1 0 .0 5 0 -0 . 0 5 -0 . 1

0

20

40

60

80

100

120

Figure 2: Largest eigen vector of R

10

2

101 10

1

101 10

1

10

1

10 10
0 100

0

10-1 10
10
0

0

10-1 10
10
-1

-1

10

-2

10-3 10
10
-2

-1

10

-2

10

-4

10

-3

0

500

10 1000 0

-3

500

10-5 10 0 10 00 0

-2

500 500

1000 1000

Figure 3: Left: Spectrum of raw data. Middle: Spectral response of FIR filter corresponding to qmax. Right: Spectral response of FIR filter when only the first 10% of qmax are used as coefficients ans spectral response of a 70th order FIR filter.

3


1. 4 1. 2 1 0. 8 0. 6 0. 4 0. 2 0 -0 . 2 -0 . 4

0

20

40

60

80

100

120

Figure 4: Result of using q noisy input data.

max

FIR filter, compared to raw data and average over 10 blocks of

1. 2 1 0. 8 0. 6 0. 4 0. 2 0 -0 . 2 -0 . 4

0

20

40

60

80

100

120

Figure 5: Result of Remez FIR filter compared with raw data and average over 10 blocks of noisy input data.

4


2 MATLAB code used for the above simulations % % % % % Explore the use of eigenanalysis techniques in post detection digital signal processing for noise rejection Potential application include spectroscopy and imaging Jon Bell 24 Jul 2000

clear % Make some fake raw data nspec = 120; ntime = 10; scale = 1.0; noisescale = 1.0; dcoff = 0.0; psi =1.0; j=1:1:nspec; k=1:1:ntime; rawdata = sin(j/3) + sin(j/5) + 1; rawdata = (rawdata.^4)/70; ncoef=12; %number of coefficient of the optimal filter to use %subplot(211) %plot(rawdata); % Turn this on for a null test with noise as an input signal %rawdata=rand(1,nspec); %Add some random noise to the data and scale it for image colourmap for i=1:ntime noisydata(i,:) = scale*(rawdata + noisescale*(rand(1,nspec)-0.5)) +dcoff; end % Make the correlation matrix and decompose it find lambdamax and qmax R = zeros(nspec); for i=1:ntime R = R + (noisydata(i,:).'*noisydata(i,:))/ntime; end %subplot(311) %plot(rawdata) %subplot(312) %plot(noisydata(1,:)) %subplot(313) %plot(mean(noisydata)) [VR,D,FLAGS] = eigs(R,1); wc=0.40; V=remez(38,[0 wc-0.04 wc+0.04 1],[1 1 0 0]); %plot(j,mean(noisydata),j,VR); %subplot(131) %S=abs(fft(rawdata,2048)); %semilogy(S(1:1000)); %subplot(132) %S=abs(fft(VR,2048)); %semilogy(S(1:1000));

5


%subplot(133) %S=abs(fft(VR(1:ncoef),2048)); %SR=abs(fft(V,2048)); %semilogy([1:1000],S(1:1000),[1:1000],SR(1:1000)); %figure(1) AR = zeros(ncoef,1); AR(1)=1; A = zeros(ntime,1); A(1)=1; for i=1:ntime cleanVR(i,:) = abs(filtfilt(VR(1:ncoef),AR,noisydata(i,:))); cleanV(i,:) = abs(filtfilt(V,A,noisydata(i,:))); end plot(j,rawdata,j,mean(noisydata),j,mean(cleanVR)) figure plot(j,rawdata,j,mean(noisydata),j,mean(cleanV))

6