算法原理参考:
https://blog.csdn.net/xvrixingkong/article/details/87860538
本仿真系统结构如下:
clear ;
close all;
clc
%% 参数设置
len = 24e2; % 每帧数据符号长度
mtype = 1; % 调制方式 1:QPSK , 2:16QAM , 3:64QAM
fame = 20*mtype; % 帧数,帧数越多误码率越精确,仿真速度越慢
M = 2*mtype ;
carrier_en = 1; %% 载波同步使能
fd = 0e3; %% 多普勒Hz
%***************************************
% 伪随机序列产生
%***************************************
% m_poly = primpoly(7,'all') ; %% 本原多项式
% 产生伪随机序列,7阶m序列,周期为127
% 利用primpoly(7,'all')指令,产生本原多项式:D^7+D^6+D^5+D^4+1,并调用编写的m_generate()函数生成伪随机序列
ff = [0 0 0 1 1 1 1]; %%低位在前
pseudoNumber = m_generate(ff);
fp = 2^length(ff) - 1; % 扩频倍数
%%
sps = 8; % % 扩频后上采样倍数
Rs = 1e5; %% 符号速率
fs = Rs*fp*sps; %% 基带采样率
% EbN0 = [4 6 8 10 12 14]; %% QPSK信噪比设置
% EbN0 = [10 12 14 16 19]; %% 16QAM信噪比设置
% EbN0 = [14 16 19 21 25]; %% 64QAM信噪比设置
EbN0 = 20;
SNR = EbN0 + 10*log10(M) - 10*log10(fs/Rs);
SNR = 20;
%% 生成bit数据
msg = randi([0 1],M*len,1); % 生成bit信号
din_data = bit2dec(msg);
din_data = [din_data din_data];
fid=fopen('qpsk_din.txt','w');
for k=1:length(din_data)
fprintf(fid,'%d \n',real(din_data(k)));
end
fclose(fid)
%% 信源编码
%% 扩频
[dsss_tx,converted,PN2] = func_dsss(msg.',pseudoNumber.',fp);
figure(1)
xx=1:length(msg(1:100));
line(xx/Rs*10e6,msg(1:100));
xlabel('时间us')
title('扩频前信息');
figure(2)
xx=1:length(dsss_tx(1:100*fp));
plot(xx/(fp*Rs)*10e6,dsss_tx(1:100*fp));
xlabel('时间us')
title('扩频后信息');
%% 扩频前频谱
[ylab,xlab] = dsss_befor_fft(msg.',mtype,Rs,4);
figure(3)
plot(xlab/1e3,ylab);
xlabel('频率kHz')
ylabel('幅度dB')
title('扩频前频谱');
%% 扩频后频谱
[ylab,xlab] = dsss_befor_fft(dsss_tx,mtype,Rs*fp,2);
figure(4)
plot(xlab/1e6,ylab);
xlabel('频率MHz')
ylabel('幅度dB')
title('扩频后频谱');
%% 扰码
sig = reshape(dsss_tx,1,[]);
[tmp,yn] = raoma(sig);
%% 调制
[tx_mod,tx_bitcoder] = modulator_new(tmp,mtype ) ;
figure(5)
plot(tx_mod(1:end),'.r');grid on;
title('调制后星座图')
figure(6)
plot(real(tx_mod),'-b');hold on;
plot(imag(tx_mod),'-r');grid on;
title('调制后时域图')
%% 调制解调验证
% [hard_bit,soft_llr,rx_bitcoder] = demodulator_soft(tx_mod.*2^15,mtype,16);
% LL = min(length(hard_bit),length(tmp));
% [E ,errBit ] = error_rate(hard_bit(1:LL),tmp(1:LL))
%% 重复多帧
head = tx_mod(1:128); %% 取N个符号作为同步帧头
tx_mod_r = repmat(tx_mod.',1,fame); %% 重复帧数
%% 2倍成型滤波
sps1 = 2; %% 上采样倍数
Alpha = 0.35; %%滚降系数
[fir_out1,fir_hn] = rcos_up(tx_mod_r,sps1,Alpha);
% figure
% fvtool(fir_out(1:1e5));
% title('成型滤波后频域图')
%% 上采样
sps2 = 4; %% 上采样倍数
[fir_out,fir_hn1] = fir_up(fir_out1,sps2);
%% IQ调制
Fs = Rs*fp*sps1*sps2; %% 50.8M
fc = Fs/4; %% 载波频率
nn=0:length(fir_out)-1;
nco = exp(1j*2*pi*fc/Fs*nn);
IQ_txout = real(fir_out).*real(nco)+ imag(fir_out).*(-imag(nco));
% IQ_txout = real(fir_out).*real(nco)+1i*imag(fir_out).*(imag(nco));
% IQ_txout = fir_out.*nco;
% figure
% fvtool(IQ_txout);
% title('IQ调制后频域图')
figure(5)
xx=1:length(tx_mod(1:1*fp));
plot(xx/(fp*Rs)*10e6,tx_mod(1:1*fp));
xlabel('时间us')
title('IQ调制前信息')
figure(6)
xx=1:length(IQ_txout(1:10*fp));
plot(xx/(fp*Rs)*10e6,IQ_txout(1:10*fp));
xlabel('时间us')
title('IQ调制后时域')
%% 加频偏
ll=1:length(IQ_txout);
nco = exp(-1j*fd/Fs*2*pi*ll);
tx_ori = IQ_txout.*nco;
%%
f1 = 10;f2=f1+1; f3=f2+1; f4=f3+1;
f5 = f4+1;f6=f5+1; f7=f6+1; f8=f7+1;
%%
for nn=1:length(SNR)
%% 加噪声
tx_agwn = awgn(tx_ori,SNR(nn));
%% IQ 解调
n=0:length(fir_out)-1;
nco = exp(-1j*2*pi*fc/Fs*n);
RX_I = real(tx_agwn).*real(nco);
RX_Q = real(tx_agwn).*(imag(nco));
rx_ori = RX_I + 1i*RX_Q;
% rx_ori = tx_agwn.*nco;
figure
fvtool(rx_ori(1:1e5));
title('IQ解制后频域图')
%% 下采样滤波抽取
rx_fir_out1= filter(fir_hn1,1, rx_ori);
rx_fir_out2 = rx_fir_out1(1:sps2:end);
%% 接收匹配滤波
rx_fir_out= filter(fir_hn,1, rx_fir_out2);
figure
fvtool(rx_fir_out(1:1e5));
title('IQ解制后频域图2')
% figure(f1-1)
% fvtool(rx_fir_out(1:1e4));
% title('接收成型滤波后频域图')
% figure(22)
% plot(rx_fir_out(1:4:end ),'.r');grid on;
% title('时间同步后星座图');
%% 时间同步
gardner_in = rx_fir_out(1:sps1/2:end);
figure (f1-1)
plot(gardner_in(1:1e4),'.r');grid on;
title('时间同步前星座图');
[gardner_out ,mu_n] = gardner_func(rx_fir_out);
figure(f1)
plot(mu_n(1:2e4));
title('时间同步收敛曲线');
%%
power = mean(abs(gardner_out(2e4:1e6)));
data_morm = gardner_out(2e4:end)./power * 0.75;
figure(f2)
plot(data_morm(end-1e5:end),'.r');grid on;
title('时间同步后星座图');
%% 载波同步
if carrier_en==1
carrier_in = data_morm;
[carrier_out, fo, c1, c2 ] = carrier_loop_new(carrier_in,fs/2,fd,mtype) ;
figure(f3)
plot(carrier_in(end-1e4:end),'.b');hold on;
plot(carrier_out(end-1e4:end),'.r');hold off;grid on;
title('载波同步前后对比')
else
carrier_out = gardner_out;
end
%% 帧同步
if mtype==3 %% 16QAM载波同步收敛慢
sync_in = carrier_out(end/2:end);
else
sync_in = carrier_out(end/4:end);
end
sync_thr = 50 ; %% 帧同步门限
[sync_out,sync_peak] = func_sync(sync_in,head,sync_thr,f4);
%%
frame_len = len*fp;
row = fix(length(sync_out)/frame_len );
for rr=1:row
phase_in = sync_out(frame_len*(rr-1)+1:frame_len*rr);
%% 相偏补偿
[phase_out,angle_head] = func_phase(phase_in,head,mtype);
figure(f5)
plot(imag(head(1:64)));hold on;
plot(imag(phase_out(1:64)),'-r');hold off;grid on;
title('帧头对比')
figure(f6)
plot(phase_in(end-1e4:end),'.b');hold on;
plot(phase_out(end-1e4:end),'.r');hold off;
title('相偏补偿前后星座图')
figure(f7)
plot(tx_mod(1:end),'.b');hold on;
plot(phase_out(1:end),'.r');hold on;grid on;
title('调制后星座图')
%% 解调
% [ReData] = demodu(mtype,phase_out);
width = 16;
standard_point_out = phase_out.*2^(width-1)*0.75/0.83;
[hard_bit,soft_llr] = demodulator_soft(standard_point_out.',mtype,width);
%% LDPC 译码
%% 解扰
[raoma_out,yn] = raoma(hard_bit);
%% 扰码前误码统计
L = min(length(sig),length(raoma_out));
[raoma_E ,raoma_errBit(rr) ] = error_rate(sig(1:L),raoma_out(1:L));
%% 解扩
rx_dsss = func_dsss2(raoma_out,pseudoNumber.',fp);
%% 误码统计
LL=M*len;
[E ,errBit(rr) ] = error_rate(rx_dsss(1:LL),msg(1:LL).');
end
nn
err_sum = sum(errBit);
E1(nn) = err_sum/(row*LL)
end
%%
figure(51)
semilogy(EbN0, E1, 'b--o');hold on;
xlabel('EbN0')
ylabel('误码率BER')
title('误码率曲线')
% figure(50)
% semilogy(xlab, y1, 'b--o');hold on;
% semilogy(xlab, y2, 'r-*');hold off;grid on;
% legend('FLOS-CMA','R-CMA')
% xlabel('样本数')
% ylabel('误码率')
如需要完整仿真程序,可私信。