Final Assignment Random Signal Processing Solution

Instructions



Write down the answers and gures in a PDF document and submit it in Brightspace. When you’re asked to give a sketch, you can make a drawing by hand, and scan it (or make a photo of it).




The report will be checked and the result of the assignment will be either PASS or FAIL. The grade for the written exam will only count in case you pass this assignment.







Signal Processing for Radio-Pulsar Naviga-tion









2.1 Introduction




A pulsar (short name for pulsating star), is a rotating neutron start that emits electromagnetic radiation. While rotating, the star emits a beam of radiation. This radiation can be observed on earth as a pulse when the beam is pointed towards earth. You can compare this with a lighthouse, where the light can only be seen when the lighthouse is pointing towards the observer.




Neutron stars have short and very regular rotational periods. As a result, the received pulses can be observed in intervals that are extremely regular. Depending on the star, the intervals (periods) between the pulses range from milliseconds to seconds. For some pulsars, the periodicity of the received pulses is more precise than an atomic clock. That makes pulsars extremely useful for many applications.




One of these (future) applications could be navigation based on pulsars, as the preciseness of the periodic pulses can be used for time-of-arrival based navigation. Among the advantages of pulsar based navigation are the fact that it is independent of satellites (like GPS) and the fact that the pulsar signals have an extremely wide frequency range, which makes it di cult to distort by potential enemies.




Besides the advantages of pulsars they carry one disadvantage, and that is the fact that the signal-to-noise ratio (SNR) of the received pulses is ex-tremely low (around -90 dB). That makes the pulsar signals very di cult to detect.




This assignment is about the detection of pulsar signals and how this can be done using concepts of random signal processing.







2















2.2 Modeling Pulsar Signals




Usually, it is assumed that the received pulsar signals are degraded by addi-tive noise that is uncorrelated with the pulsar signal. That is,




y(t) = x(t) + n(t);
(1)



where t indicates the (continuous) time index, x(t) the pulsar signal, n(t) the noise and y(t) the received noisy pulsar signal. The pulsar signal is assumed to be deterministic, while the noise n(t) is assumed to be a stochastic sta-tionary and ergodic uncorrelated signal with zero mean, i.e., E[N(t)] = 0, and variance N2 . For many pulsars, the period of the pulsar signal x(t), say T0, is known by approximation. This information is obtained by analyzing extremely long recordings of received noisy pulsar signals. By performing speci c signal processing operations on these recordings it is also possible to estimate one period of a pulsar signal x(t). A method that is traditionally often used is called epoch folding. With epoch folding one divides the noisy signal in frames of length T0 and averages the data over these periods. As the noise is assumed to be uncorrelated, it will cancel out when enough averaging operations are performed.




Exercise 1




Sketch the autocorrelation function of the noise process N(t)










2.3 Epoch folding




The Epoch folded signal is given by




1
K 1






Xk


ze(t) = K
(2)
y(t + kT0);




=0





where K denotes the number of frames over which the signal is averaged. Substituting Eq. (1) into (2) we obtain




1
K 1




1
K 1












X










Xk






ze(t) = K




K




(3)
x(t + kT0) +
n(t + kT0) ;






k=0










=0










|




{z




}
|




{z




}








x0(t)




n(t) with n(t) ! 0
as K!1





3

























































Figure 1: Example of epoch folding.




where x0(t) indicates one complete period of the pulsar signal and n(t) the averaged noise. As the noise is an uncorrelated zero-mean process, n(t) will average out if K is su ciently large. An example of this procedure is given in Fig. 1. To get an impression how well the pulse can be detected in the noise, one can compute the SNR. Let the SNR of the measured noisy signal, computed over one period, be given by






0


R
T=00 x0(t)2dt
1




SN R = 10 log10
@
E
t
T0
n(t)2dt
iA
:
(4)






hR
t=0









The SNR of the epoch folded signal is given by






0


R
T0 x0(t)2dt
1




SN Re = 10 log10
@
E


A
:
(5)
t
T0
n(t)2dt










=0












hRt=0


i




Exercise 2




Compute how much the SNR of the epoch folded signal has increased com-pared to the original SNR if K = 10 averaging operations are performed (that is, compute SN Re SN R).




Exercise 3




How large should K be, in order to obtain an SNR increase of 90 dB using epoch folding? Assume that one period of the pulsar signal takes 0.5 seconds. How many hours of data do we then need to increase the SNR with 90 dB?










4















From the blackboard website you can download a Matlab data- le named pulsar data.mat. Download this le and load it using the command







load pulsar data,







This data le contains a synthetically generated example pulsar signal in the variable x, the noise in variable n, the mixed signal in variable y and a variable called template that we will use later on. These signals are sampled at 10 kHz.




Exercise 4

Determine the periodicity T0 in samples by observing signal x. Also compute the SNR of the original signal using signals x and n.




Exercise 5

Write an m- le [ze]=epoch folding(signal, T0, K);




This m- le takes a signal as input, together with the period T0 and the num-ber of frames K over which we average. The output is an epoch folded signal, i.e., one period of the signal.




Exercise 6

Use the m- le [ze]=epoch folding(signal, T0, K); to perform epoch folding in order to detect the pulse in the pulsar signal. Use for input vari-able signal the noisy signal y and use for input variable T0 the period that you determined in exercise 4. For K, use the values that should give you an SNR increase of 10, 20, and 30 dB. Give the corresponding values for K and plot the three resulting epoch folded signals for the three di erent SNRs. Can you detect the position of the pulse in the signal?













2.4 Matched ltering




The epoch folding approach estimates the signal of interest. However, in order to have a su ciently high output SNR such that detection is easy, K should be very large for these extremely low SNRs. Moreover, we are mainly interested in the position of the pulse in the signal, and not so much in the signal itself. Furthermore, for many pulsars, the signal x is known as researchers performed epoch folding techniques on extremely long data-




5















records. This means that templates of one period of the signal x(t), i.e., x0(t), are available. The availability of such a template simpli es detection of the pulse a lot as additional information on the signal is used. The data- le that you loaded also contains one period of a template in the variable template.




A well-known method that makes use of such a template is the matched-ltering approach. The matched lter has a very speci c impulse response h(t). Under the signal model y = x + n we can write the output of the matched lter by the following convolution




zm(t) =
Z +1
y(s)h(t
s)ds


(6)
Z +1


=
x(s)h(t
s)ds + Z +1
n(s)h(t s)ds
(7)
=
zx(t) + zn(t);




(8)



where zx(t) and zn(t) are the signals that result when the impulse response of the matched lter is applied to signal x(t) and n(t), respectively. To derive the impulse response of the matched lter, the per sample ratio between signal power and noise variance is maximized, that is,




( ) =
jzx( )j2
:
(9)
E [jzn( )j2]












Using the Cauchy-Schwartz inequality, it can be shown that ( ) is maxi-mized by choosing h(t) as h(t) = cx( t + ), that is, a mirrored and shifted version of x(t) (due to t and in the argument of x( t + ), respectively) and where c is an arbitrary constant. This impulse response can also be seen as a template of the signal of interest. By convolving the noisy signal y with this impulse response, a correlation between the noisy signal y and the known template of x is essentially computed, that is,




+1
zm(t) = c y(s)x( t + + s)ds: (10)




For now we consider the case for which = 0, that is,




+1
zm(t) = c y(s)x( t + s)ds: (11)










6















2.5 Relation Between Matched Filtering and Epoch Folding




By limiting the integration boundary of Eq. (11) to an integer multiple of the pulse period T0, say KT0, it is possible to show that the matched ltering can be implemented as the correlation of one pulse period of the template and the estimate ze(t) obtained by epoch folding, that is,

KT0
zm(t) =


y(s)x(
t + s)ds


(12)
=
0


K k=0
y(s + kT0)!
x( t + s)ds
(13)
Z0




T0
1
K 1












X









T0
=
ze(s)x( t + s)ds;
(14)


0





where constant c is set to c = K1 From this we can conclude that the matched lter can indeed be implemented as the correlation between one pulse period of the template, i.e., x( t + s), and the estimate ze(t) that is obtained by epoch folding.







By substituting Eq. (3) into Eq. (14) this can also be written in terms of x0(t) and n(t),

ZT0 ZT0

zm(t) = x0(s)x( t + s)ds + n(s)x( t + s)ds: (15)

0 0




Exercise 7

Assume that x0(t) consists of a perfect pulse as shown in Fig. 2. Make a sketch of R0T0 x0(s)x( t + s)ds. Explain the shape of the sketch.




The matched lter output zm(t) given by Eq. (15) still contains some noise that is expressed by the second term. The expected value of the statistical process that underlies zm(t), i.e., E [Zm(t)], can be shown to be




E [Zm(t)] = T0RX0X ( t):
(16)



Exercise 8




Derive the expression in Eq. (16) and explain which statistical assumptions you make in order to derive this.










7
























x0(t)



















1







0.5







0

0 0.2 0.4 0.6 0.8 1




t



Figure 2: Example of one period of x(t).







In Matab we work with data that is sampled at a certain sampling fre-quency fs. To implement Eq. (15) in Matlab, the integrals in Eq. (15) need to be written as summations and instead of continuous time we use sampled time. This leads to




1
N 1






X


zm(l) =


ze(m)x( l + m)
(17)
fs




m=0


1
=
T0
X
ze(m)x(
l + m)




(18)
















N
m=0
























=
T0
N 1
x0(m)x(
l + m) +
T0
N 1
n(m)x( l + m); (19)




X






X




N
m=0




N
m=0

















where the relation 1 = T0 is used.




fs N

Exercise 9




For large N, i.e., N ! 1, zm(l) converges to




zm(l)
T0RX0X(
l) + T0RNX ( l)
(20)
=
T0RX0X(
l):
(21)



Which property related to estimation of autocorrelation functions should the underlying processes have such that the approximation from Eq. (18) to (20) is valid.




Eq. (18) resembles a correlation estimator between the epoch folded signal ze(l) and the template x( l). Implementing this estimator can be done












8















by writing Eq. (18) in terms of a matrix-vector multiplication zm = xze, i.e.,

2
zm(1)
3


T0
2
x( 1)
x(0)
: : :
x(N
2)32
ze(1)
3


zm(0)








x(0)
x(1)
: : :
x(N
1)
ze(0)


6
...
7
=


6
...


. . .
...
76
...
7
N
N +2)
6zm(N
1)7




6
x( N + 1) x(
: : :
x(0)
76ze(N
1)7
6


7




6








76


7
4


5




4








54


5
As x(m) is periodic with N samples, we can replace x(m) by x(m mod N). Leading to zm = xmod N ze, that is,

2
zm(1)
3


T0
2
x(N
1) x(0)
: : :
x(N
2)
32
ze(1)
3
6
zm(0)
7




6
x(0)
x(1)
: : :
x(N
1)
76
ze(0)
7(22):
...
=


...


. . .
...


...
N
x(2)


6zm(N
1)7




6
x(1)
: : :
x(0)


76
ze(N
1)7
6


7




6










76


7
4


5




4










54


5
To construct the matrix xmod N you can make use of the Matlab command toeplitz.




Exercise 10

Write an m- le [zm]=matched filter[ze,x], which performs the above matrix-vector multiplication in order to compute the output of the matched lter. Make use of the Matlab command toeplitz to implement above ma-trix structure. To test the implementation use ze = [x0(0); :::; x0(N 1)]T and compute zm = xze. Make a plot of the result. Is this in line with your sketch from Exercise 7?







Exercise 11




We will now use this implementation of the matched lter to determine the pulse location of the pulsar signal. First generate epoch folded signals using the function epoch folding on the noisy signal y for the K-values that you found in Exercise 6. Use each of these epoch folded signals as input ze in the function matched filter. For each K-value, make a plot of the epoch-folded signal and the matched lter output. Which value of K do you need such that you can detect the pulse using epoch folding? Which value of K do you need to detect the pulse if you combine epoch-folding and the matched lter? Give the relative position of the detected pulse in terms of seconds compared to the template.
















9