数字通信——MSK调制解调,维特比解码下误码率与信噪比之间的关系
标签: matlab 数字通信 算法
前言
一个实现数字通信下msk调制与解调的MATLAB仿真,以及利用维特比算法实现维特比解码。通过维特比解码来实现最大似然解码的一些个人见解,最后得出信噪比与误码率之间的图像。
一、MSK是什么?
在这里插入图片描述
二、仿真代码
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
2.维特比算法代码
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
图像
信噪比与误码率的函数图像:
智能推荐
jQuery基础总结(三)
1.1 复习jQuery操作DOM 1.2 元素操作 1.2.1 高度和宽度 1.2.2 坐标值 $(“div”).offset(); // 获取或设置坐标值 设置值后变成相对定位 $(“div”).position(); // 获取坐标值 子绝父相 只能读取不能设置 1.2.3 滚动条(滚动事件) $(“div”).scro...
配置STP定时器
配置STP定时器 1、 PC1 ping PC2 2、配置STP定时器 在4台交换机上配置使用STP,并配置S1为该二层网络中根交换机,SW2为备份根交换机。 配置完成查看定时器的默认值】 在PC4上一直ping PC2,持续发送ICMP报文,并进行连通性测试。 在SW1上修改STP的Forwad Delay时间为2000cs,默认为1500CS,cs代表百分之一秒,只有在根交换机上配置才会生效。...
美图网上太多难以筛选,教你用Python挑选最合适的
前几天,極光同学写了篇下载王者荣耀皮肤的文章,可以轻松的获取各种英雄背景图,甚是激动,也想将桌面背景换成漂亮的,不过我对王者荣耀不感冒(曝露年龄啦),时常会被必应搜索主页的背景图所震撼,于是想到从必应获取桌面壁纸,并且排除掉自己不喜欢的图片,应该是个不错的主意,说干就干 问题分析 必应每天都会有新的壁纸,大都是自然风光、人文地理等等,非常漂亮,在页面上点击右键,保存背景图片,就能简单的保存下来。但...
MyBatis入门简单操作
MyBatis入门简单操作 1.Mybatis简介 MyBatis 是一款优秀的持久层框架,它支持定制化 SQL、存储过程以及高级映射。MyBatis 避免了几乎所有的 JDBC 代码和手动设置参数以及获取结果集。MyBatis 可以使用简单的 XML 或注解来配置和映射原生信息,将接口和 Java 的 POJOs(Plain Ordinary Java Object,普通的 Java对象)映射成...
基础C语言知识串串香9☞C语言复杂表达式
文章参考微信公众号[嵌入式软件学习圈] 四、C语言复杂表达式 4.1、在表达式中,要看符号的优先级和结合性。 4.2、在理解内存时,内存0地址在最底下,至上地址逐渐增加。 4.3、int *p;是定义的一指针变量p,而int (*p)[4];也是一个指针变量p;也可以这样想:凡是遇到(*p)什么的判断他是指针后,就可以说他是指针变量,包括函数指针。 4.4、一个函数int max(int a, i...
猜你喜欢
github上传本地项目
github学习:如何从本地把项目上传到github&&如何把github项目通过clone复制下来,详细教程 一、第一步---注册一个Github账号 首先要在GitHub上创建一个帐号,可以去官方网站注册一个账号。 前提:本地安装一个git 本人github:https://github.com/saucxs 二、git安装 下载地址:http://m...
GPON介绍及华为OLT网关注册配置流程
一、GPON介绍 1.GPON简介 随着技术的发展,光纤变得“便宜又好用”,因此FTTx(FiberTo The X,光纤接入)作为新一代宽带解决方案被广泛应用,它为用户提供高带宽、全业务的接入平台。同时,FTTH(FiberTo The Home,光纤到户)更是被称为是最理想的业务透明网络,是接入网发展的最终方式。 那么,FTTx是如何实现的呢?在多种方案中,PON(Pa...
spring 5.x 学习笔记
目录 第一章 引言 1. EJB存在的问题 2. 什么是Spring? 3. 设计模式 4. 工厂设计模式 总结 第一章 引言 1. EJB存在的问题 .运行环境苛刻 .代码移植差 总结:EJB是重量级框架 2. 什么是S...
leetcode+华为笔试题-java实现返回一个整数数组中最大子数组的和
方法一:暴力枚举 定义一个最大值max初始化一个很小的数,定义一个变量sum表示求和值,遍历数组元素,从第一个元素开始,依次相加,如果和sum比最大值max大就将sum赋值给最大值。然后再来一个循环控制从第i个数组元素开始求和,直到n. 时间复杂度:O(n^2) 方式二:贪心法 因为每次求和都是将i前面的元素相加,会出现重复的情况,如果出现相加完之后是负数,说明这时我就要继续遍历找到第一个正数,将...
SpringCloud G版 eureka-consumer-ribbon 服务消费 ribbon
eureka-consumer-ribbon 系列目录 一、父工程 二、服务注册 三、服务发现 四、服务消费 五、服务消费 ribbon 六、服务消费 openfeign 1.建module 2.勾依赖 3.启动类@EnableDiscoveryClient 4.改配置 5.Controller 6.启动测试 7.访问http://localhost:9002/hello 8.ribbon完成 下...