¹⁰² views
ubuntu2004

Kernel: Octave

In [0]:

Lab session V : Identification of signals of known shape using correlation

Introduction

In this lab, we will implement a detector able to to discriminate a signal of known shape either buried in noise or mixed among others over a given temporal window. The detector relies on the concept of cross-correlation (also referred to as a “sliding scalar product”) between signals, a method which allows for the detection or the recognition of a signal with a known shape.

The name "sliding scalar product" comes from the fact that for a given t, the expression of the cross-correlation $R_{xy}(t)=\int x(\tau) y(\tau - t) d\tau = \int x(\tau +t) y(\tau) d\tau$ between two real signals is merely the scalar product of $y(t)$ with $x(\tau +t)$ , a replica of $x(\tau)$ shifted by t (to the left).

Widely used, the detection by means of correlation is at the basis of useful technical applications like radars and sonars, which are able to cope with noisy environments and weak reflected signals. Correlation also allows sending several signals simultaneously on the same channel, a process called multiplexing, and to decipher them at the reception as well.

Problem 5.1 : Implementing the running correlation algorithm

Implement a function returning as a vector the result of the numerical cross-correlation $R_{xy}[k] = \sum x[n+k] y[n]$ between two input vectors x and y, y[n] being the reference signal (i.e. the - short - signal of known shape), x[n] being the signal within which the occurrence(s) of signal y are to be detected.

The running correlation process may be thought of as a sampled system which provides a sample of the correlation signal r[n] for every input sample x[n]. Hence, the operation may be performed “in real time”, in the sense that every new input sample (coming from a radar for instance) yields a new correlation sample. This process can be done sequentially and thus could be made very quickly (“running correlation algorithm”), depending on the processing power.

Load the contents of the file Correl_data.mat into your Cocalc environment ( load Correl_data ). Display in a graph the signals x[n] and y[n] which are contained within the loaded file. Both of them are sampled at Te=0.01 s. Identify the positions (i.e. the indexes) of the three occurrences of the rectangular impulsion y[n] mixed with the shot noise within x[n].
Using your function, calculate the cross-correlation between x[n] and y[n] and check whether your function runs properly or not.
What happens if the order of the arguments x and y is inverted when calling the function ?

In [3]:

function R = mycorr(x,y,Ts)

nx = length(x);
ny = length(y);

N = nx + ny

buffer = zeros(1,ny);

x = [x,buffer];
buffer_index = 2:length(buffer);
R = zeros(1,N);

for i = 1:N
    buffer = [buffer(buffer_index),x(i)];
    R(i) = sum( buffer .* conj(y)) * Ts;
end



end

In [7]:

load correl_data

In [6]:

%exc 1
Ts = 0.01; 



Nx = length(x);
tx = 0:Ts:(Nx*Ts-Ts);

Ny = length(y);
ty = 0:Ts:(Ny*Ts-Ts);

subplot(2,1,1)
grid on
plot(tx,x)
grid on 
subplot(2,1,2)
grid on 
plot(ty,y)


%the rect pulse can be identified at 30,50 and 80 seconds

Out[6]:

In [7]:

R = mycorr(x,y,0.01);
plot(abs(R))
grid on

Out[7]:

N = 11001

In [8]:

%exc 3 
%reversing the order of arguments reverses the time of R
R = mycorr(y,x,0.01);
plot(abs(R))
grid on

Out[8]:

N = 11001

In [5]:

%method to find maxima in the correlation
function out = find_max(signal,num) % num is number of maxima that should be returned 

out = [];


for i = 1:num
   
    maximum = max(signal);
    n = find(signal == maximum);
    out = [out,n];
    signal(n) = 0;
end

end

In [10]:


R = mycorr(x,y,0.01);
maxima_guess = zeros(3,2);

maxima_guess(1,:)=[3500,4000];
maxima_guess(2,:)=[5000,6000];
maxima_guess(3,:)=[8300,10000];
index_maxima=[];

for i = maxima_guess'

    a = i(1);
    b = i(2);
    index_maxima =[index_maxima,(find_max(abs(R(a:b)),1)+a-1)];
    
end

%exact times of occurance
tx(index_maxima)-ty(end)

Out[10]:

N = 11001
ans =

   29   49   79

In [11]:

 [maximum,index] = max(R)

Out[11]:

maximum = 2.1435
index = 3901

In [12]:

%exact occurences
plot(R)
grid on 
hold on

for i = 1:length(index_maxima)

    tmp=index_maxima(i);
    plot([tmp,tmp],[-3,3],'LineWidth',1,'r');
    hold on

end

Out[12]:

Problem 5.2 : Detection of radar impulsions

The file correl_data.mat contains two other signals: radar_pulse (the reference radar impulsion) and radar_received, which stands for the signal detected by the radar system and gathers one or several reflections of the reference impulsion mixed with a strong white noise.

Draw both signals on the same graph (rescale radar_pulse if it is too small). Calculate the energy of the radar impulsion. Can you identify the occurrences of the reflected impulsions within the detected signal ?
Use the function you have conceived above to compute the cross-correlation $R_{xy}(t)$ between these signals and draw the result in a new graph. Where are the reflected impulsions located ?
In a real radar system, the detection criterion is given by the comparison between the correlated signal and an alert threshold. We first set the alert threshold equal to half the energy of the reference impulsion. Chose a realistic threshold and calculate a binary alarm signal whose samples are defined as: $R_{xy}(t) > Threshold$ . Append this signal (in red) to the former graph.
The sampling frequency of the radar signal is $10^7$ samples/s. Deduce the distance of each object behind each reflected impulsion.
If the detected signal is too noisy, detection errors may occur (wrong alerts). To assess the errors rate of the detector, compute the cross-correlation between radar_pulse and the signal radar_noise, using the same detection threshold as before. radar_noise mimics a signal without reflected impulsions. Draw the result of the cross-correlation.
How many times does the cross-correlation exceed the alert threshold ( hint : the instruction: val=sum(corr>threshold)) ? Assess the error rate in terms of events/s. How many times is the cross-correlation smaller than the threshold minus the energy of the reference impulsion ? Assess the non-detected event rate in terms of events/s.
Consider questions 5 and 6 again with two other values of the alert threshold, the first lower, the second larger than the one used so far. Comment your results.

In [3]:

%exc 1
R = mycorr(radar_pulse,radar_pulse,1/10^7);
Energy = R(round(length(R)/2))

Out[3]:

N = 198
Energy = 3.7500e-06

In [15]:

length(support)

Out[15]:

ans = 199

In [16]:

length(R)/2
support = -98:99;

Out[16]:

ans = 99

In [17]:

plot(support,R)

Out[17]:

In [15]:

plot(radar_received)

Out[15]:

In [8]:

%exc 1
Ts=1/10^7
t=0:Ts:(length(radar_received)*Ts-Ts);
nx = length(radar_received)
ny = length(radar_pulse)
nzeros = round((nx-ny)/2)-1

zero = zeros(1,nzeros);
pulse2 = [zero,radar_pulse,zero,[0]];


plot(t,radar_received)
hold on
plot(t,pulse2*5,'r')
grid on 
legend('received Signal','amplified radar pulse')
%the signal can not be indentified with the eye

Out[8]:

Ts = 1.0000e-07
nx = 5000
ny = 99
nzeros = 2450

In [11]:

%exc 2
R = mycorr(radar_received,radar_pulse,Ts);

Out[11]:

N = 5099

In [12]:

%find maxima
maxima_guess(1,:)=[500,900];
maxima_guess(2,:)=[900,1100];
maxima_guess(3,:)=[1500,2000];
maxima_guess(4,:)=[4000,5000];
index_maxima=[];

for i = maxima_guess'

    a = i(1);
    b = i(2);
    index_maxima =[index_maxima,(find_max(abs(R(a:b)),1)+a-1)];
    
end

In [13]:

plot(R)
grid on 
hold on 

for i = 1:length(index_maxima)

    tmp=index_maxima(i);
    plot([tmp,tmp],[-10^(-5),10^(-5)],'LineWidth',1,'r');
    hold on

end

ylim([-10^(-5),10^(-5)])

legend('correaltion','position of maxima')

Out[13]:

In [21]:

%exc 3
threshold = Energy/2

tmp = find( R > threshold );

binary = zeros(length(R));
binary(tmp) = threshold; % random value just to make sure both signals fit easily in the same plot 
a = binary(1:5099); %for some reason if i dont do this extra step the plot of binary always causes cocalc to crash

plot(R)
hold on 
plot(a)
grid on 
legend('Corrrelation','Binary signal')
%from this plot its not that clear to me if the actual threshold value is supposed to be smaller or bigger than E/2
%I assume a higher threshold value is required since there are 4 distinct peaks and with the value E/2 there is one additional peak

Out[21]:

threshold = 1.8750e-06

In [42]:

threshold_2 = 2e-6

tmp = find( R > threshold_2 );

binary = zeros(length(R));
binary(tmp) = threshold_2;
a=binary(1:5099);

plot(R)
hold on 
plot(a)
grid on 
legend('Corrrelation','Binary signal 2')

Out[42]:

threshold_2 = 2.0000e-06

In [14]:

%exc4 
%We use the maxima found in the correlation before to find the time values of the occurrences in the received signal 
maxima_index_signal = index_maxima-length(radar_pulse)
Ts = 1/10^7
T = maxima_index_signal * Ts % time delay until the signal is received 
v_light = 300000000
distances = T*v_light/2 %divided by 2 since the signal is traveling the distance twice

Out[14]:

maxima_index_signal =

    566    933   1820   4191

Ts = 1.0000e-07
T =

   5.6600e-05   9.3300e-05   1.8200e-04   4.1910e-04

v_light = 3.0000e+08
distances =

    8490   13995   27300   62865

In [44]:

%exc5
R_noise = mycorr(radar_pulse,radar_noise,Ts);
t5=0:Ts:Ts*length(R_noise)-Ts;
plot(t5,R_noise)
hold on
plot(t5,t5*0+threshold_2)
grid on 
%This shows that the threshold value of 2e-6 is not high enough becaue the correlation with the noise exceeds this value several times
%based on this 3e-6 is a good estimation 
val = sum(R_noise>threshold_2)
errorr_rate=val/length(R) % approx 1.9% of the samples are detected wrongly 


val2 = sum(R_noise<threshold_2-Energy)
error_rate = val2/length(R)
legend('correlation pulse and white noise','threshold value')

Out[44]:

N = 5099
val = 97
errorr_rate = 0.019023
val2 = 160
error_rate = 0.031379

In [45]:

subplot(2,1,1)
t6=0:Ts:Ts*length(R)-Ts;
threshold_3 = 1e-6

tmp = find( R > threshold_3 );

binary = zeros(length(R));
binary(tmp) = threshold_3;
a=binary(1:5099);

plot(t6,R)
hold on 
plot(t6,a)
grid on 
legend('Corrrelation','Binary signal')
title(['Threshold =',num2str(threshold_3)])


subplot(2,1,2)
threshold_4 = 4e-6

tmp = find( R > threshold_4 );

binary = zeros(length(R));
binary(tmp) = threshold_4;
a=binary(1:5099);

plot(t6,R)
hold on 
plot(t6,a)
grid on 
legend('Corrrelation','Binary signal')
title(['Threshold =',num2str(threshold_4)])

Out[45]:

threshold_3 = 1.0000e-06
threshold_4 = 4.0000e-06

Problem 5.3 : Multiplexing digital signals using CDMA

Introduction

In certain standards of mobile cell phones, multiple users transmit digital (binary) signals across a common transmission channel, a small portion of the electromagnetic spectrum in the microwave range. Several signals may be transmitted at any time, completely uncoordinated, and thus be present at the same time on the channel. A process called multiplexing allows the proper decoding of the message bits transmitted by each user.

A correlation-based detector may be used for this purpose, provided that each user pair use a common and specific digital binary code (a private key), i.e. a short sequence of 1 and 0, embodied for transmission in an analog sequence of rectangular impulses (e.g. a small rectangular burts at +5V to code 1 and at -5V to code 0).

This code or its inverse is used for transmission of the binary signal : to send a sequence of message bits, the user concatenates these positive (coding a 1) and negative (coding a 0) versions of the code signal into a sequence of code signals which are sent to the communication channel. Other users transmit their own message bits in the same fashion, except that, of course, they use different private keys.

A communication system of this form is called a code-division, multiple-access (CDMA) system or a direct-sequence, spread-spectrum (DSSS) system.

To decode the message, the receiver correlates the multiplexed signal with the private key. Evaluated after each code duration (bit-time), a running correlation will evaluate the resemblance or anti-resemblance of the "noisy" signal (where several user massages may superpose !) with the private key. The actual guess of the transmitted sequence is merely based on the sign of the correlator ouput. The uncertainty on the received values is higher when the absolute value of the correlator output is close to zero.

Simulation of CDMA detection

The file correl_data also contains the variable dsss (direct-sequence spread-spectrum) which mimics a multiplexed signal built as the superposition of four different binary signals, sent by four user pairs (one user pair means a sender and a receiver).

We want to determine the signal transmitted through a given user pair, so called bit-sequence, in terms of a sequence of bits (0 or 1). Every bit among a user pair’s bit-sequence is coded by a code sequence : bit 1 stands for the user’s code sequence itself and bit 0 for its opposite. User pairs code sequences are all different and may also be of different lengths. All code sequences are assumed to be properly synchronized at the beginning of dsss. Here, the codes sequences of the first two user pairs cs1 and cs2 are assumed to be known.

Draw the graphs of dsss, cs1 and cs2 within the same figure. Calculate the energy of cs1 and cs2.
Use the function you have conceived at exercise 1 to compute the cross-correlation between dsss and the longest code sequence, cs1. Call the resulting vector cor1, and draw its graph in a new figure.
The decoding of the first user’s bit-sequence must be performed by extracting the proper samples of cor1. Indeed, one only needs the samples of the cross-correlation standing for the temporal coincidence between dsss and cs1 (in other words, the only relevant information is the vector of the scalar products between every entire bit-sequence and the user pair code sequence). Superpose these samples to the former graph. To make them more visible, draw them using the instruction stem and cor1 with plot.
The samples extracted from the cross-correlation can be used to derive the transmitted user’s bit-sequence. If their value is positive, the received bit is assumed to be “1”, or “0” otherwise. Following this, decode the bit-sequence transmitted by user 1 (Hint : it consists of ten bits, the first three being 0 1 1).
Repeat this procedure to decode the bit-sequence transmitted by user 2 (Hint : it consists of 17 bits). The code sequence cs2 being shorter than cs1, its energy is smaller. Hence, the bit error probability at the reception is large. Identify the most doubtful bits.

In [4]:

%exc1
subplot(3,1,1)
plot(dsss)
subplot(3,1,2)
plot(cs1)
subplot(3,1,3)
plot(cs2)
R = mycorr(cs1,cs1,1);
Energy1 = R(round(length(R)/2))
R = mycorr(cs2,cs2,1);
Energy2 = R(round(length(R)/2))

Out[4]:

N = 200
Energy1 = 100
N = 120
Energy2 = 60

In [5]:

1 > 4

Out[5]:

ans = 0

In [16]:

%exc2-4 
cor1 = mycorr(dsss,cs1,1);
plot(cor1)


sig1_peaks = zeros(1,length(cor1)); 
index = length(cs1):length(cs1):length(cor1)-length(cs1);%index of values where the signal is shifted an integer time of its length 
sig1_peaks(index) = 1;%set these values to one
bit_seq_1 = cor1(index)>0%if the correlation is bigger than 0 at this points its a 1 if not a 0


hold on 
stem(sig1_peaks.*cor1)
grid on

Out[16]:

N = 1100
bit_seq_1 =

  0  1  1  1  1  1  0  0  0  0

In [23]:

%exc5
cor2 = mycorr(dsss,cs2,1);
plot(cor2)


sig2_peaks = zeros(1,length(cor2)); 
index = length(cs2):length(cs2):length(cor2);%index of values where the signal is shifted an integer time of its length 
sig2_peaks(index) = 1;%set these values to one
bit_seq_2 = cor2(index)>0%if the correlation is bigger than 0 at this points its a 1 if not a 0


hold on 
stem(sig2_peaks.*cor2)
grid on

Out[23]:

N = 1060
bit_seq_2 =

  0  1  1  0  1  0  0  0  0  0  1  0  1  1  1  1  0

In [25]:

%most doubtful bits 

stem(1./abs(cor2(index)),'r')
grid on 

%the 17th bit is the most doubtful one 
%3,7 and 11 are the next must doubtful ones

Out[25]:

Conclusions

Complete this cell with your personal conclusions about what you've learned. Give an estimation on the time you've spent (outside the alloted time for lab sessions) realizing this lab session.

In [0]:

Lab session V : Identification of signals of known shape using correlation

Introduction

Problem 5.1 : Implementing the running correlation algorithm

Problem 5.2 : Detection of radar impulsions

Problem 5.3 : Multiplexing digital signals using CDMA

Introduction

Simulation of CDMA detection

Conclusions

Product

Resources

Company