1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
clear all
close all
clc

%-------------------------------------------------------------------------%
%                                PART 1.1                                 %
%-------------------------------------------------------------------------%

N = 1e4; % 10^5
A = 1;
T = 0.01;

% 生成一个符号序列{-1,1}
b = (sign(randn(1,2*N))+1)/2;
for i=1:1:length(b) 
    if b(i)==0
        b(i)=-1;
    end
end

% 应用递归方程对QPSK中的符号进行变换
% symbols
xQ_1 = -1;
x_1 = 1;
xI_n = zeros(1,N);
xQ_n = zeros(1,N);
for n=1:1:N
    if n==1
        xI_n(n) = -xQ_1*x_1;
        xQ_n(n) = -xI_n(n)*b(n);
    elseif n~=1 
        xI_n(n) = -xQ_n(n-1)*b(2*(n-1));
        xQ_n(n) = -xI_n(n)*b(2*n-1);
    end
end

xI_n; %  同相分量
xQ_n; %  正交分量
z_n = xI_n + 1i*xQ_n; % Add in - phase and quadrature component to create z_n

% 计算误码率
iter = 1e2; 
SNR_dB = 5;
SNR_lin = 10.^(SNR_dB/10); % 10log_10(SNR_lin)
BER_appr = 0;

for p=1:1:iter
    n = randn(1,N) + 1i*randn(1,N);
    y_n = A*T*z_n + sqrt(A^2*T^2/SNR_lin)*n;
    y_Re = sign(real(y_n)); % 当出现负数时返回-1，否则返回1
    y_Im = sign(imag(y_n));
    BER_appr = BER_appr + sum(y_Im ~= imag(z_n)) + sum(y_Re ~= real(z_n));
end
BER = BER_appr/(N*iter);

D = ['BER = ', num2str(BER)];
disp('------------------------------------------------------------------------------------------------------------------------------------');
disp('An approximation for BER - Bit Error Rate with SNR = 5 dB is:');
disp(D);
disp('------------------------------------------------------------------------------------------------------------------------------------');

SNR_dB_B = 5:12; % Values of SNR in dB
SNR_lin_B = 10.^([5:12]/10); % Values of SNR in decimal

BER_B = zeros(1,length(SNR_lin_B));
BER_appr_B = zeros(1,length(SNR_lin_B));

for k=1:1:length(SNR_lin_B)
    for p=1:1:iter
        n = randn(1,N) + 1i*randn(1,N);
        y_n = A*T*z_n + sqrt(A^2*T^2/SNR_lin_B(k))*n;
        y_Re = sign(real(y_n)); % returns -1 when a negative number occur and 1 otherwise
        y_Im = sign(imag(y_n));
        BER_appr_B(k) = BER_appr_B(k) + sum(y_Im ~= imag(z_n)) + sum(y_Re ~= real(z_n));
    end
    BER_B(k) = BER_appr_B(k)/(2*N*iter);
    Theory_BER(k) = 0.5*erfc(sqrt(SNR_lin_B(k))); % theoretical BER
end

disp(['SNR in dB : ' num2str(SNR_dB_B)]);
disp(['BER Measured : ' num2str(BER_B)]);
disp(['BER Theoretical : ' num2str(Theory_BER)]);
disp('------------------------------------------------------------------------------------------------------------------------------------');

figure(1)
semilogy(SNR_dB_B,BER_B,'b-*') 
xlabel('SNR (dB)')
ylabel('BER - Bit Error Rate Measured')
title('BER - Bit Error Rate vs SNR_{dB}')
grid on

figure(2)
semilogy(SNR_dB_B,BER_B,'b-*') 
hold on
semilogy(SNR_dB_B,Theory_BER,'m-o')
xlabel('SNR (dB)')
ylabel('BER - Bit Error Rate Measured & Theoretical')
title('BER - Bit Error Rate vs SNR_{dB}')
legend('BER - Measured', 'BER - Theoretical')
grid on

%-------------------------------------------------------------------------%
%                                PART 1.2                                 %
%-------------------------------------------------------------------------%

SNR_dB_C = 5;
SNR_lin_C = 10.^(SNR_dB_C/10);

% Creating the symbols x(n) in set {-1,1}
b = (sign(randn(1,N))+1)/2;
for i=1:1:length(b) 
    if b(i)==0
        b(i)=-1;
    end
end

% Constructing the vectors s_1 & s1
% s_1 is for s^(-1)
% s1 is for s^1
s_1 = [-(2*A*sqrt(T)*1i)/pi; (A*sqrt(T)*sqrt(pi^2 - 4))/pi];
s1  = [A*sqrt(T); 0];

beta = (A^2*T)/SNR_lin_C; % From equation (16)
var = 2*beta;
BER_Viterbi = 0; %initialising BER

for p=1:1:iter
    n1_n = sqrt(var)*(randn(1,N) + 1j*randn(1,N)); % Generating the values of the noise vector n1_n
    n2_n = sqrt(var)*(randn(1,N) + 1j*randn(1,N));
    phase(1) = 0;
    for n=1:1:N
        phase(n+1) = phase(n) + b(n)*pi/2;
        if b(n) == -1 
            r_n(:,n) = s_1.*exp(1i*phase(n)) + [n1_n(n); n2_n(n)]; % if x(n) == -1
        elseif b(n) == 1
            r_n(:,n) = s1.*exp(1i*phase(n)) + [n1_n(n); n2_n(n)];  % if x(n) == 1
        end
    end
    x_opt = Viterbi(N,s1,s_1,r_n);
    BER_Viterbi = BER_Viterbi + sum(b~=x_opt); 
end
BER_Viterbi = BER_Viterbi/(2*N*iter);

disp('An approximation for BER - Bit Error Rate with Viterbi aglgorithm and SNR = 5 dB is:');
disp(['BER = ', num2str(BER_Viterbi)])
% disp(['BER Theoretical : ' num2str(Theory_BER)]);
disp('------------------------------------------------------------------------------------------------------------------------------------');

% Creating the BER diagram for SNR = 5,6,7,8,9,10,11,12
SNR_dB_E = 5:12;
SNR_lin_E = 10.^(SNR_dB_E/10);
BER_Viterbi_E = zeros(1,length(SNR_lin_E));
BER_Viterbi_appr_E = zeros(1,length(SNR_lin_E));

% Constructing the vectors s_1 & s1
s_1 = [-(2*A*sqrt(T)*1i)/pi; (A*sqrt(T)*sqrt(pi^2 - 4))/pi];
s1  = [A*sqrt(T); 0];

for k=1:1:length(SNR_lin_E)
    for p=1:1:iter
        n1_n = sqrt(A^2*T/SNR_lin_E(k))*(randn(1,N) + 1j*randn(1,N)); % Generating the values of the noise vector n1_n
        n2_n = sqrt(A^2*T/SNR_lin_E(k))*(randn(1,N) + 1j*randn(1,N));
        phase(1) = 0;
        for n=1:1:N
            phase(n+1) = phase(n) + b(n)*pi/2;
            if b(n) == -1 
                r_n(:,n) = s_1.*exp(1i*phase(n)) + [n1_n(n); n2_n(n)]; % if x(n) == -1
            elseif b(n) == 1
                r_n(:,n) = s1.*exp(1i*phase(n)) + [n1_n(n); n2_n(n)];  % if x(n) == 1
            end
        end
        x_opt = Viterbi(N,s1,s_1,r_n);
        BER_Viterbi_appr_E(k) = BER_Viterbi_appr_E(k) + sum(b~=x_opt);
        Theory_BER_Viterbi(k) = 0.5*erfc(sqrt(SNR_lin_E(k))); % theoretical BER
    end
    BER_Viterbi_E(k) = BER_Viterbi_appr_E(k)/(2*N*iter);
end

disp(['SNR in dB : ' num2str(SNR_dB_E)]);
disp(['BER Measured : ' num2str(BER_Viterbi_E)]);
disp(['BER Theoretical : ' num2str(Theory_BER_Viterbi)]);
disp('------------------------------------------------------------------------------------------------------------------------------------');


figure(3)
semilogy(SNR_dB_E,BER_Viterbi_E,'b-*') 
xlabel('SNR (dB)')
ylabel('BER - Bit Error Rate Measured ')
title('BER - Bit Error Rate vs SNR_{dB} - Viterbi Algorithm')
grid on

figure(4)
semilogy(SNR_dB_E,BER_Viterbi_E,'b-*'),
hold on
semilogy(SNR_dB_B,BER_B,'g-*'),
hold on
semilogy(SNR_dB_B,Theory_BER,'m-o')
xlabel('SNR (dB)')
ylabel('BER - Bit Error Rate Measured & Theoretical')
title('BER - Bit Error Rate vs SNR_{dB}')
legend('BER - Viterbi','BER - Measured', 'BER - Theoretical','Location','southwest')
grid on

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
function [x_opt] = Viterbi(N, s1, s_1, r_n)

%-------------------------------------------------------------------------%
%                                 FORWARD                                 %
%-------------------------------------------------------------------------%

pointer_pi(1) = -1;     
pointer_0(1) = -1;

for n=1:1:N
    if n==1
        w_3pi2(n) = real(r_n(:,n)'*s_1*exp(1i*0)); % 第一步
        w_pi2(n) = real(r_n(:,n)'*s1*exp(1i*0));   
    elseif n~=1
        if mod(n,2)==0 % Even bits
            % Even symbols can end up ONLY with phase pi or 0
            
            % From 3pi/2 to pi with symbol -1
            % From 3pi/2 to 0 with symbol +1
            w3pi2_pi(n) = real(r_n(:,n)'*s_1*exp(1i*3*pi/2));
            w3pi2_0(n) = real(r_n(:,n)'*s1*exp(1i*3*pi/2));
            
            % From pi/2 to pi with symbol +1
            % From pi/2 to 0 with symbol -1
            wpi2_pi(n) = real(r_n(:,n)'*s1*exp(1i*pi/2));
            wpi2_0(n) = real(r_n(:,n)'*s_1*exp(1i*pi/2));

            % The cost may be the weight(3pi/2, pi) + the weight of the
            % last symbol 3pi/2 due to the memory property of the phase.
            % The cost may be the weight(pi/2, pi) + the weight of the
            % last symbol pi/2 due to the memory property of the phase.
            total_cost_1 = w3pi2_pi(n) + w_3pi2(n-1);
            total_cost_2 = wpi2_pi(n) + w_pi2(n-1);
            [w_pi(n),pointer_pi(n)] = max([total_cost_1 0 total_cost_2 0]);
            
            % The cost may be the weight(3pi/2, 0) + the weight of the
            % last symbol 0 due to the memory property of the phase.
            % The cost may be the weight(pi/2, 0) + the weight of the
            % last symbol 0 due to the memory property of the phase.
            total_cost_1 = w3pi2_0(n) + w_3pi2(n-1);
            total_cost_2 = wpi2_0(n) + w_pi2(n-1);
            [w_0(n),pointer_0(n)] = max([total_cost_1 0 total_cost_2 0]);
            
        elseif mod(n,2)~=0 % Odd bits
            % Odd symbols can end up ONLY with phase 3pi/2 or pi/2
            
            % From 0 to 3pi/2 with symbol -1
            % From 0 to pi/2 with symbol +1
            w0_3pi2(n) = real(r_n(:,n)'*s_1);
            w0_pi2(n) = real(r_n(:,n)'*s1);
            
            % From pi to 3pi/2 with symbol +1
            % From pi to pi/2 with symbol -1
            wpi_3pi2(n) = real(r_n(:,n)'*s1*exp(1j*pi));
            wpi_pi2(n) = real(r_n(:,n)'*s_1*exp(1j*pi));
            
            % The cost may be the weight(pi, 3pi/2) + the weight of the
            % last symbol pi due to the memory property of the phase.
            % The cost may be the weight(0, 3pi/2) + the weight of the
            % last symbol 0 due to the memory property of the phase.
            total_cost_1 = wpi_3pi2(n) + w_pi(n-1);
            total_cost_2 = w0_3pi2(n) + w_0(n-1);
            [w_3pi2(n),pointer_3pi2(n)] = max([0 total_cost_1 0 total_cost_2]);

            % The cost may be the weight(pi, pi/2) + the weight of the
            % last symbol pi due to the memory property of the phase.
            % The cost may be the weight(0, pi/2) + the weight of the
            % last symbol 0 due to the memory property of the phase.
            total_cost_1 = wpi_pi2(n) + w_pi(n-1);
            total_cost_2 = w0_pi2(n) + w_0(n-1);
            [w_pi2(n),pointer_pi2(n)] = max([0 total_cost_1 0 total_cost_2]);
        end
    end
end

%-------------------------------------------------------------------------%
%                                 BACKWARD                                %
%-------------------------------------------------------------------------%

% 只保留给出最大权重和的路径
if mod(n,2)~=0
    [~,route(N+1)] = max([w_3pi2(n) 0 w_pi2(n) 0]);
elseif mod(n,2)==0
    [~,route(N+1)] = max([0 w_pi(n) 0 w_0(n)]);
end

route(1) = 4;
for n=N:-1:1
    if n~=1
        if mod(n,2)==0 % Even symbols
            [~,p] = max([0 w_pi(n) 0 w_0(n)]);
            enter = [0 pointer_pi(n) 0 pointer_0(n)];
            route(n) = enter(p);
        elseif mod(n,2)~=0 % Odd symbols
            [~,p] = max([w_3pi2(n) 0 w_pi2(n) 0]);
            enter = [pointer_3pi2(n) 0 pointer_pi2(n) 0];
            route(n) = enter(p);
        end
    end
    route;
    
    % Restoring the symbols 
    if route(n)-route(n+1)==-1
        x_opt(n) = -1;
    elseif route(n)-route(n+1)==1
        x_opt(n) = 1;
    elseif route(n)-route(n+1)==-3
        x_opt(n) = 1;
    elseif route(n)-route(n+1)==3
        x_opt(n) = -1;
    end  
    
end

end