BPSK-DSSS收发系统

MATLAB仿真

目录

• MATLAB仿真
  • 目录
  • 系统实现
    • 系统框图
  • 仿真代码
    • 生成m序列
    • 扩频码与映射序列相乘
    • 根升余弦滤波
    • 上变频（乘以载波）
    • 物理信道（添加噪声）
    • 下变频（乘以载波）
    • 解扩+匹配滤波
    • 采样+判决
    • 画出系统输入输出序列对比
  • 误码率曲线仿真
    • 代码
    • 输出图像
  • 附录代码
    • 代码1
    • 代码2

系统实现

系统框图

唯一与BPSK不同的是，对极信号与扩频码进行相乘。

仿真代码

生成m序列

% 生成m序列
data_len = 10;               %数据长度 
n = 3;                       %阶次
p = 2^n - 1;                 %循环周期
ms = idinput(p,'prbs');
mg = repmat(ms',1,data_len); %重复m序列

figure;
subplot(121);
stairs(ms);title('m序列时域波形');
ylim([-1.5 1.5]);
c = xcorr(ms,'normalized');
subplot(122);
plot(c);title('m序列自相关函数');

扩频码与映射序列相乘

• 方法一：（方便理解）

  data_len = 10;                       %数据长度 
  data = randi([0 1],1,data_len);      %生成0/1序列
  data_bio = (data-1/2)*2;             %映射，0变-1
  s = rectpulse(data_bio,p);           %序列重复p倍
  y = s.*mg;                           %信号与扩频码相乘

  figure;
  subplot(211);
  stairs(s);title('扩频之前的序列波形');
  ylim([-1.5 1.5]);
  subplot(212);
  stairs(y);title('扩频之后的序列波形');
  ylim([-1.5 1.5]);

  

• 方法二：（矩阵重组）

  思路，生成的m序列为列向量ms，乘以之前的映射序列data_bio；生成一个数组，然后利用reshape函数重组数组，得到y（相乘之后的序列）；

  reshape函数按列取元素

  n = 3;          %阶次
  p = 2^n - 1;     %循环周期
  ms = idinput(p,'prbs');
  m_1 = ms.*data_bio;
  y1 = reshape(m_1,[1,p*data_len]);
  figure;
  subplot(211);
  stairs(s);title('扩频之前的序列波形');
  ylim([-1.5 1.5]);
  subplot(212);
  stairs(y1);title('扩频之后的序列波形');
  ylim([-1.5 1.5]);

  

根升余弦滤波

osr = 8;                      % 8倍上采样
y_up = upsample(y,osr); % 上采样
fir_rcos = rcosdesign(alfs,6,osr,'sqrt');
data_bx = conv(y_up,fir_rcos);            %根升余弦滤波

上变频（乘以载波）

fs = 16000;                   % 采样率 
fc = 2000;                    % 载波频率
carrier = cos(2*pi*fc/fs*(0:length(data_bx)-1)); %载波
s_up = carrier.*data_bx;                    %上变频（乘以载波）

物理信道（添加噪声）

r_n = awgn(s_up,10,'measured');               %添加噪声

下变频（乘以载波）

r_down = r_n .* carrier;                    %下变频

解扩+匹配滤波

mg_up = upsample(fliplr(ms'),osr);        %先把序列进行反序
fir_m = conv(mg_up,fir_rcos);             %本地广义成形

r_bx = conv(r_down,fir_m);                %接收匹配滤波器

采样+判决

recv_data = (r_bx(97:56:end-104)) > 0;
%recv_data = (r_bx(osr*6+p*osr-osr+1:p*osr:end-osr*6-p*osr)) > 0;

本代码中滤波器的长度*fir_rcos**为49，7位m序列，8倍上采样，即56个点；所以总延时为49-1+56=104，由于最后面是7个0，所以卷积处理时不会体现，即最开始为104-7=97，当然也可以取序列长度为1，寻找第一个最佳相关峰；一般地，如下：recv_data = (r_bx(osr*6+p*osr-osr+1 : p*osr : end-osr**6-p*osr)) > 0;

画出系统输入输出序列对比

figure;
subplot(211);
stem(data,'filled');title('原始数据');
subplot(212);
stem(recv_data,'filled');title('输出序列');

由于序列长度较少，所以理论上出错的可能性较小，即输出序列应该和输入序列完全一致，即可认为系统设计初步验证无误。

误码率曲线仿真

代码

% bpsk-dsss 误码率仿真（输出bpsk及bpsk-dsss）
clc;clear all;close all;
fs = 16000;                   % 采样率 
fc = 2000;                    % 载波频率
%symbol_rate = 1000;           %符号速率
% Rb = 1000;                    % bit速率
alfs = 0.35;                   % 滚降系数
data_len = 100000;               %数据长度
data = randi([0 1],1,data_len);      %生成0/1序列
data_bio = (data-1/2)*2;             %映射，0变-1

%% bpsk-dsss相关操作
% 生成m序列
n = 3;          %阶次
p = 2^n - 1;     %循环周期
ms = idinput(p,'prbs');
mg = repmat(ms',1,data_len);  %重复m序列

s = rectpulse(data_bio,p);
y = s.*mg;        %信号与扩频码相乘

osr = 8;                      % 8倍上采样
s_up_d = upsample(y,osr); % 上采样
fir_rcos = rcosdesign(alfs,6,osr,'sqrt');
data_bx_d = conv(s_up_d,fir_rcos);            %根升余弦滤波

carrier_d = cos(2*pi*fc/fs*(0:length(data_bx_d)-1)); %载波
s_carrier_d = carrier_d.*data_bx_d;                    %上变频（乘以载波）

%% bpsk(无扩频操作)
s_up = upsample(data_bio,osr);   %直接上采样
data_bx = conv(s_up,fir_rcos);
carrier = cos(2*pi*fc/fs*(0:length(data_bx)-1)); %载波
s_carrier = carrier.*data_bx;

% 求误码率
errbit = zeros(1,10);
errbit_d =zeros(1,10);
snr_d = zeros(1,10);
snr = zeros(1,10);
i=1;
for ebn0 = 1:10
  snr(i) = ebn0 + 10*log10(2) - 10*log10(osr); %SNR和EbN0的转化
  %snr_d(i) = ebn0 + 10*log10(2) - 10*log10(osr) - 10*log10(p);
    snr_d(i) = snr(i) - 10*log10(p);
  %物理信道（加性高斯白噪声）
  r_noise = awgn(s_carrier,snr(i),'measured');
  r_noise_d = awgn(s_carrier_d,snr_d(i),'measured');
  
  %下变频（乘以载波）
  r_signal_down = r_noise.*carrier;
  r_signal_down_d = r_noise_d.*carrier_d;
  
  %匹配滤波器 根升余弦滤波
  r_bx = conv(r_signal_down,fir_rcos);
  %解扩加匹配滤波
  mg_upsample = upsample(fliplr(ms'),osr) ;
  fir_m = conv(mg_upsample,fir_rcos);              %本地广义成形
  r_bx_d = conv(r_signal_down_d,fir_m);                %接收匹配滤波器
  %判决
  recv_data = (r_bx(osr*6+1:osr:end-osr*6)) > 0;
  errbit(i) = sum(abs(recv_data-data))/data_len;
  recv_data_d = (r_bx_d(97:56:end-96)) > 0;
  errbit_d(i) = sum(abs(recv_data_d-data))/data_len;
  i = i + 1;
end

ebn0 = 1:10;
ber_theory = 0.5*erfc(sqrt(power(10,ebn0/10)));%理论误码率
figure
semilogy(ebn0,ber_theory,'-or');
grid on;hold on;

semilogy(ebn0,errbit,'-sb');
grid on;hold on;
title('BPSK误码率曲线图');

semilogy(ebn0,errbit_d,'-xg');
grid on;hold on;
xlabel('EbN0(dB)');ylabel('BER')
title('BPSK-DSSS误码率曲线图');
legend('理论误码率曲线','bpsk仿真误码率曲线','bpsk-dsss仿真误码率');

输出图像

附录代码

代码1

clc; clear; close all
%:本程序用于实现QPSK直接扩频系统
%========================基本参数===========================%
Rs = 100e6; %符号速率
sps = 4; %上采样倍数
fs = sps * Rs; %采样率
alfs = 0.35; % 滚降系数
data_len = 1e3; %数据长度
M = 4; %调制阶数
offset = pi / 4; %偏移量
rng(9);
EbN0 = 0:2:10;
snr = EbN0 - 10 * log10(sps) - 10 * log10(2 ^ n) + 10 * log10(2);
% snr = 100;
ber = zeros(1, length(snr));
for i = 1:length(snr)
    bit_error = 0;
    k = 0;
    while (bit_error > 100 && k * data_len > 1e4)
        %========================调制==============================%
        bit_send = randi([0 1], 1, data_len); %bit数据
        data = bit2int(bit_send.', log2(M)).';
        % data = randi([0 M - 1], 1, data_len); %随机数据
        data_mod = pskmod(data, M, offset); %PSK调制
        % scatterplot(data_mod); title('QPSK-Constellation Diagram');
        %========================扩频==============================%
        n = 10; %m序列级数
        C = 435; %八进制反馈系数
        tap_set = [1 0 0 0 0 0 0 1 0 0 1]; %二进制
        m_sequence = 1 - 2 * [m_serial_generator(tap_set) 1]; %生成扩频序列
        data_matrix = (m_sequence.') .* data_mod; %扩频矩阵
        data_dsss = reshape(data_matrix, 1, []);
        %=======================成形滤波器==========================%
        Fir_rcos = rcosdesign(alfs, 8, sps, 'sqrt');
        data_dsss_up = upsample(data_dsss, sps);
        data_dsss_shape = conv(data_dsss_up, Fir_rcos);
        %=========================噪声=============================%
        data_awgn = awgn(data_dsss_shape, snr(i), 'measured');
        %=======================解扩+匹配===========================%
        m_sequence_up = upsample(fliplr(m_sequence), sps);
        Fir_match = conv(m_sequence_up, Fir_rcos); %本地广义匹配滤波器
        data_matched = conv(data_awgn, Fir_match); %解扩+匹配
        %========================解调===============================%
        data_matched_down = data_matched(4125:sps * 1024:end - 4124);
        data_demod = pskdemod(data_matched_down, M, offset);
        bit_recv = int2bit(data_demod.', log2(M)).';
        bit_error = bit_error + biterr(bit_send, bit_recv);
        disp(['Bit Error = ', num2str(bit_error)]);
        k = k + 1;
    end
    ber(i) = bit_error / (k * data_len);
    disp(['BER = ', num2str(ber(i))]);
end
%==========================画图=============================%
ber_theory = berawgn(EbN0.', "psk", M, "nondiff");
figure
semilogy(EbN0, ber, '-*g'); grid on; hold on;
semilogy(EbN0, ber_theory, '-r'); grid on; hold on;

%===================扩频码生成函数===========================%
function [m] = m_serial_generator(tap_set)
    L = 2 ^ (length(tap_set) - 1) - 1;
    x = [zeros(1, (length(tap_set) - 2)) 1];
    for i = 1:L
        m(i) = x(end);
        for j = 1:length(tap_set) - 1
            sum_vector(j) = tap_set(j + 1) * x(j); %#ok<*AGROW>
        end
        sum_x = mod(sum(sum_vector), 2);
        x = [sum_x x(1:end - 1)];
    end
end

代码2

clc; clear; close all
%:本程序用于实现QPSK直接扩频系统
%========================基本参数===========================%
Rs = 100e6; %符号速率
sps = 4; %上采样倍数
fs = sps * Rs; %采样率
alfs = 0.35; % 滚降系数
data_len = 1e4; %数据长度
M = 4; %调制阶数
offset = pi / 4; %偏移量
rng(9);
%========================调制==============================%
bit_send = randi([0 1], 1, data_len); %bit数据
% bit_send = repmat([0 1],[1 data_len/2]);
data = bit2int(bit_send.', log2(M)).';
% data = randi([0 M - 1], 1, data_len); %随机数据
data_mod = pskmod(data, M, offset); %PSK调制
scatterplot(data_mod); title('QPSK-Constellation Diagram');
%========================扩频==============================%
n = 10; %m序列级数
C = 435; %八进制反馈系数
tap_set = [1 0 0 0 0 0 0 1 0 0 1]; %二进制
% m_sequence = [idinput(2^n-1,'prbs').' 1];
m_sequence = 1 - 2 * [m_serial_generator(tap_set) 0]; %生成扩频序列
data_matrix = (m_sequence.') .* data_mod; %扩频矩阵
data_dsss = reshape(data_matrix, 1, []);
% data_mod_repmat = rectpulse(data_mod,1024);
% mg = repmat(m_sequence,[1,length(data_mod)]);
% data_dsss = data_mod_repmat.*mg;
%=======================成形滤波器==========================%
Fir_rcos = rcosdesign(alfs, 8, sps, 'sqrt');
data_qpsk = conv(upsample(data_mod, sps), Fir_rcos);
data_dsss_up = upsample(data_dsss, sps);
data_dsss_shape = conv(data_dsss_up, Fir_rcos);
%=========================噪声=============================%
snr = 10;
data_qpsk_awgn = awgn(data_qpsk, snr, 'measured');
data_awgn = awgn(data_dsss_shape, snr, 'measured');
%=======================解扩+匹配===========================%
data_qpsk_matched = conv(data_qpsk_awgn, Fir_rcos); %匹配滤波
m_sequence_up = upsample(fliplr(m_sequence), sps);
Fir_match = conv(m_sequence_up, Fir_rcos); %本地广义匹配滤波器
data_matched = conv(data_awgn, Fir_match); %解扩+匹配
%========================解调===============================%
data_qpsk_matched_down = data_qpsk_matched(33:sps:end - 32);
scatterplot(data_qpsk_matched_down); title('SNR=10dB QPSK无扩频');
data_matched_down = data_matched(4125:sps * 1024:end - 4124);
scatterplot(data_matched_down); title('SNR=10dB QPSK有扩频');
data_qpsk_demod = pskdemod(data_qpsk_matched_down, M, offset);
bit_qpsk = int2bit(data_qpsk_demod.', log2(M)).';
bit_error = biterr(bit_send, bit_qpsk);
disp(['QPSK Bit Error = ', num2str(bit_error)]);
ber = bit_error / data_len;
disp(['QPSK BER = ', num2str(ber)]);
data_demod = pskdemod(data_matched_down, M, offset);
bit_recv = int2bit(data_demod.', log2(M)).';
bit_error = biterr(bit_send, bit_recv);
disp(['QPSK-DSSS Bit Error = ', num2str(bit_error)]);
ber = bit_error / data_len;
disp(['QPSK-DSSS BER = ', num2str(ber)]);
%========================频谱分析===========================%
fft_n = 16384; %FFT点数
f_1 = (-fft_n / 2:fft_n / 2 - 1) ./ fs;
f_2 = (-fft_n / 2:fft_n / 2 - 1) * 1024 ./ fs;
data_qpsk_awgn_fft = fft(data_qpsk_awgn, fft_n);
data_awgn_fft = fft(data_awgn, fft_n);
figure;
plot(f_1, 20 * log10(abs(fftshift(data_qpsk_awgn_fft)) / fft_n), 'r'); hold on;
plot(f_2, 20 * log10(abs(fftshift(data_awgn_fft)) / fft_n), 'b'); grid on;
xlabel('Frequency (Hz)'); ylabel('Power (dB)'); legend('QPSK', 'QPSK-DSSS');
%===================扩频码生成函数===========================%
function [m] = m_serial_generator(tap_set)
    L = 2 ^ (length(tap_set) - 1) - 1;
    x = [zeros(1, (length(tap_set) - 2)) 1];
    for i = 1:L
        m(i) = x(end);
        for j = 1:length(tap_set) - 1
            sum_vector(j) = tap_set(j + 1) * x(j); %#ok<*AGROW>
        end
        sum_x = mod(sum(sum_vector), 2);
        x = [sum_x x(1:end - 1)];
    end
end