BER result in MATLAB

Question

I build a pam-2 modulation, then I make a pulse shaping with half sine (matched filter). Then I send it through the AWGN channel. At the end I do down sampling and demodulation. But i have problem with plotting BER. I don't understand what I'm doing wrong:

clc;
clear;
N=1e4;
N2 = 1e2; 
M = 2;
range = 0:10;
error =  zeros(1,length(range)); %BER

%
% half sine
Rc = 1e3;       % Chip rate
T = 1/Rc;      % inverse of chip rate
Tc = 0.5* T;
Fs = 2e3;      % sampling frequency
dt = 1/Fs;
over = Fs/Rc;   % sampling factor
sps = 10;  

time = 0:dt/sps:2*T;
half_Sine = sin(pi*time/(2*T)).^3; 

%% BER

for i = 1:length(range)
    for n = 1:N2 
        
% Modulation
        x=randi([0 M-1],N,1);
        h_mod = pammod(x,M);       
        over_data=upsample(h_mod,over); 
        txSig = conv(over_data,half_Sine, 'same'); 
% AWGN  
        Ps = mean(((txSig)).^2);
        Sigma = sqrt(Ps  * 10^(-range(i)/10) / 2);
        Noise = randn(length(txSig), 1) * Sigma;
        rx_SIG = Noise + txSig; 

% Downsample
     down = rx_SIG(1:over:end);
     
%  Demodulation 
        hDemod  =  pamdemod(down,M);
% Errors
       error(i) = error(i)+...
           sum(hDemod~=x) / length(hDemod);

    end
    
     BER = error/n;
end
figure(1);
grid on
semilogy(range,BER);
title('BER');

Update:I need to build a ber from 10^1 to 10^6

Answer 1

1.- Your script up and running :

There are many code lines turned comments that I used to get to run the start script. I have left them along, that should be used while debugging only.

close all;clear all;clc;

N=1024;
M = 2;

% pulse : half sine
Rc = 1e3; % [b/s] chip rate
T = 1/Rc; % [s/b] inverse of chip rate
Tc = 0.5* T;
Fs = 2e3; % [Hz] sampling frequency
dt = 1/Fs; % [s]  
ov1 = Fs/Rc;   % sampling factor
sps = 10;  
dt2=dt/sps

% single pulse time reference
t = [0:dt2:2*T];  % [s]

% signals usually have a known preable
% a heading known sequence that helps receivers
% extract with relative ease when pulse peaks take place
% the sync preamble has to be long enough to acquire
% the sampling interval.
% For simplicity here I just make sure that the first bit of the signal is always 1 
nsync=64 % length sync header
x=[ones(1,nsync) randi([0 M-1],1,N-nsync)];  % signal : data

xm=reshape(x([nsync+1:end]),[(N-nsync)/8 8])
L1=sum(repmat(2.^[7:-1:0],size(xm,1),1).*xm,2); % string to check received bytes against, to measure BER

over_data = pammod(x,M); % signal : PAM symbols

% err1 =  zeros(1,length(x)); %BER

% single pulse
pulse_half_Sine = sin(pi*t/(2*T)).^3; 
figure;plot(t,pulse_half_Sine)
grid on;xlabel('t');title('single pulse')

% A=[.0001:.0001:1];
A=[1e-3 5e-3 1e-2 5e-2 .1 .5 1 10 100];
% A=[5e-2:1e-2:2];
rng1 = [1:numel(A)];  % amount S power levels to check for BER

% to use power on requires a reference impedance
% usually assumed 1, in accademic literature, but then when attempting to correlated
% BER with used Watts if R0=1 the comparison is not correct.
R0=50 % [Ohm]


%% BER measuring loop
% k=1

% Logging BER for different signal  power levels,
% Logging signal power levels, it was useful when getting script up and running
BER=[]; 
SNR_log=[]; 

figure(1)
ax1=gca

for k = 1:length(rng1)

% generating signal
x2=2*(x-.5); % [0 1] to [-1 1]
 S=[]
 for k2=1:1:numel(x)
     S=[S  A(k)*x2(k2)*pulse_half_Sine];
 end
 
Ps = mean(S.^2)/R0; % signal power
        
% adding AWGN 
% you are making the noise proportional to the signal
% this causes BER not to improve when signal power up
% sigma1 = sqrt(Ps  * 10^(-rng1(k)/10) / 2); % not used
sigma1=.1
noise1 = randn(length(S), 1) * sigma1;

Pn=mean(noise1.^2)/R0;
        
rx_S = noise1' + S; % noise + S

% Downsample
% this downsampling is an attempt to sync received signal
% to the time stamps where pulse peaks are expected
% but it does not work
%      dwn1 = rx_SIG(1:ov1:end);
     
%  Demodulation 
% because the sampling times of the previous line are not
% centered pamdemod doesn't work either
%         hDemod  =  pamdemod(dwn1,M);

% the missing key step : conv on reception with expected pulse shape
rx2_S=conv(pulse_half_Sine,rx_S);

rx2_sync=conv(pulse_half_Sine,rx_S([1:1:nsync*numel(pulse_half_Sine)]));

% removing leading samples that only correspond to the 
% pulse used to correlate over received signal
rx2_S([1:numel(pulse_half_Sine)-1])=[];

% syncing 
[pks,locs]=findpeaks(abs(rx2_sync),'NPeaks',nsync,'MinPeakHeight',A(k)/2);
% [pks,locs]=findpeaks(abs(rx2_S));

% x3(find(pks<.1))=[]; % do not use sign results close to zero
% locs(find(pks<.1))=[];

%  t0=dt2*[0:1:numel(rx2_sync)-1];
%  figure;plot(t0,rx2_sync);hold on
%  plot(t0(locs),pks,'bo')

% 5 header pulses needed, 1st and last header samples are null, not needed
n01=find(pks<.2*A(k));
if ~isempty(n01) peaks(n01)=[];locs(n01)=[]; end
%  pks([1 end])=[];locs([1 end])=[];
%  plot(t0(locs),pks,'rs')

% since we know there have to be 5 leading pulses to be all ones
% we extract the sampling interval from this header

nT2=round(mean(diff(locs)))

% t3=dt2*[0:1:numel(rx2_S)-1];
% figure;plot(t3,abs(rx2_S));
% hold on
% xlabel('t')
% plot(t3(locs),pks,'ro')

% nt3=[1:nT2:numel(x)*nT2];  % sampling times
% plot(t3(nt3),max(pks)*ones(1,numel(nt3)),'sg')

x3=sign(rx2_S([1:nT2:N*nT2])); % only N bits expected so only sample N times

% x3 [-1 1] back to [0 1] otherwise when comparing v3 against x
% a roughtly 50% of bits are always going to be wrong, which comes
% from the signal statistics
x3=.5*(1+x3);

% sampling
% x3=sign(rx_S(locs));

% making sure x3 and x same length
% x3(find(pks<.1))=[]; % do not use sign results close to zero

SNR=Ps/Pn;
SNR_log=[SNR_log Ps];

x3_8=reshape(x3([nsync+1:end]),[(N-nsync)/8 8])
Lrx=sum(repmat(2.^[7:-1:0],size(x3_8,1),1).*x3_8,2);

err1 = sum(L1~=Lrx) / length(Lrx);
BER = [BER err1];

end

%% BER(S)
figure(1);
plot(rng1,BER);
grid on
title('BER/rng1');

This is not BIT ERROR RATIO as we know it and as it is used in all sort of quality measurements.

Note that Bit Error Ratio is not the same as Bit Error Rate despite both terms commonly used alike.

A rate implies and amount/seconds a velocity, speed.

BER as used commonly used to measure signal quality is a RATIO, not a rate.

BER = correct_bits/total_bits , but it's not as simple as this, as I am going to show.

For instance note that worst BER obtained with your script with a quick fix doesn't reach above 0.5 (!?) BER certainly reaches 1 when message not 'getting-there'.

I believe the following points are important for you to understand how BER really works.

2.- BER was completely flat for really dispare signal power levels

In an earlier working script not shown even using pulse amplitude A=100, and low noise mean(noise1)=-7.36e-04 about 1/3 of the received symbols are erroneous while figure;plot(rx_S) shows a rather clean signal, no riding ripple, no sudden changes ..

The 1/3 errorenous bit were not corrupted by channel noise but it was already in the transmitted signal. I have spaced each pulse enough to avoid overlapped pulses.

Adjacent pulses need at least 2ms to avoid overlapping.

This is without considering doppler.

Heavily overlapping symbols is what happens when command conv is used on a train of pulses generated the way you did :

S = conv(over_data,A(k)*pulse_half_Sine, 'same');

3.- You started with 1e4 data bits treated as 1e4 modulation symbols

But your transmitted-received time signal also showed length 1e4 time samples, cannot be, way too few time samples.

The time reference of over_data and pulse_half_Sine should not be the same. Nyquist; signal is currupted beyond recovery if only 2 samples er cycle of let's say carrier modulating pulses.

I tried

h_mod = pammod(x,M);       
over_data=upsample(h_mod,ov1); 
S = conv(h_mod,A(k)*pulse_half_Sine, 'same'); % modulated signal

h_mod = pammod(x,M);       
S = conv(h_mod,A(k)*pulse_half_Sine, 'same'); % modulated signal

S = conv(over_data,A(k)*pulse_half_Sine, 'same'); % modulated signal

and none of these 3 got the expected BER showing whether the signal is strong or weak.

4.- It turns out command upsample is for discrete-time models

sys = tf(0.75,[1 10 2],2.25)
L = 14;
sys1 = upsample(sys,L)

not to directly interpolate a signal to, for instance, double the amount of samples as it seems you attempted.

5.- This is how the transmitted signal (before noise added) should look like

t2=dt2*[0:1:numel(S)-1];
figure;plot(t2,S);
grid on;xlabel('t');title('transmitted signal before noise')

t3=dt2*[0:1:numel(rx2_S)-1];
[pks,locs]=findpeaks(abs(rx2_S))
figure;plot(t3,rx2_S);
hold on
xlabel('t')
plot(t3(locs),pks,'ro')

6.- The chosen pulse is not particularly strong against AWGN

The main reason being because it's a baseband pulse. not modulated, and on top of this only has positive values.

Convolution efficiency highly improves when modulating the pulse, the positive and negative pulse samples to be found across each pulse increases robustness when attempting to decide whether there's pulse or just noise.

For instance Chirp pulses are a lot stronger.

7.- To measure BER : Use bytes, constellation points, coded symbols, but not bare bits

Measuring BER with bare bits, or more broadly speaking, using a random test signal with fixed statistical moments BER is constrained to whatever mean and var assinged to signal and/or mean var from noise in absence of with weak signal.

Rewording, testing for BER with bare bits counting, when weak or no signal BER is actually measuring the noise the signal was trying to avoid.

Roughly 50% of the received bits, regardless of signal or noise, the way you are attempting BER measurement, will always hit what are apparently correct bits : false positives.

To avoid these false positives following I show how to measure BER against expected caracters.

N=1024
..
nsync=64 % length sync header
x=[ones(1,nsync) randi([0 M-1],1,N-nsync)];  % signal : data

Now x is 1024 and the initial 64 bits are for syncing only, leaving N-sync for message.

Let's check BER against let's say L1 the expected sequence of bytes

xm=reshape(x([nsync+1:end]),[(N-nsync)/8 8])
L1=sum(repmat(2.^[7:-1:0],size(xm,1),1).*xm,2);

L1 is checked against Lrx generated with x3_8 the message part of x3 the demodulated symbols

8.- the upsampling downsampling didn't work

this downsampling on reception

dwn1 = rx_SIG(1:ov1:end);

was an attempt to sync received signal to the time stamps where pulse peaks are expected but it didn't not work.

Because the sampling times were not centered pamdemod didn't work either.

9.- Use sync header to calculate sampling interval

I only convolve the nsync (64) initial bits

rx2_sync=conv(pulse_half_Sine,rx_S([1:1:nsync*numel(pulse_half_Sine)]));

These pulses allow a reliable calculation of nT2 the sampling interval to check along the rest of the received frame.

I obtain nT2 with

[pks,locs]=findpeaks(abs(rx2_sync),'NPeaks',nsync,'MinPeakHeight',A(k)/2);

There's need for further conditioning but basically locs already has the necessary information to obtain nT2 .

10.- This is the graph obtained

when no signal BER = 1 and when signal strength high enough PAM signals show good `BER' ending to 0.

When refining A step, by this meaning making it smaller, one gets the following

BER testers are often plugged to base stations upon setup and left a few hours or even days recording, and such testers do not record bare bit errors, bytes, constellation points, and even frames are checked.

11.- BER/SNR BER/EbN0 not against signal only

BER is usually plotted against SNR (analog signals) or Eb/N0 (digital signals) not just against signal amplitude or signal power.

12.- The Communications Toolbox is an add-on

This toolbox adds the following support functions: pammod pamdemod genqammod genqamdemod, yes pammod and pamdemod use genqammod genqamdemod respectively.

These functions are not available unless the Communications Toolbox is installed.

for BER simulations try Simulink, there are already available BER examples.