信噪比为36dB，进行的操作。
1 %%%%%%%%%%%%%%%%%%%%% 下行同步时频估计 %%%%%%%%% 2 %%%%%%%%%%%%%%%%%%%%% sync.m %%%%%%%%% 3 %%%%% data:2022年2月6日 author:wupingping %%%%%%%%%% 4 5 %%%%%程序说明 6 %%先生成接收信号以及同步序列，完成了BPSK调制,得到经成型滤波处理后的接收信号， 7 %%对接收信号进行加频率偏移和时延，通过信道之后，利用本地同步序列进行时频偏估计 8 %%可以通过频率估计均方误差以及定时均方误差来判断是否合理???不确定 9 10 %通信*具体内容如下： 11 %调制方式：BPSK 编码方式：无 12 %滚降因子：0.35 13 %噪声：线性高斯白噪声 14 %本地同步序列：经BPSK调制过的m序列 15 16 %%%%%仿真环境 17 %软件版本：MATLAB R2021b 18 19 20 clc; 21 clear all; 22 close all; 23 24 %********************** 程序主体 ************% 25 26 %%%%%%%%%%%%%%%%% 参数设定 %%%%%%%%%%%%%% 27 bit_rate = 12000; %比特率 28 % symbol_rate = 12000; %符号率 29 sps = 12; %每个符号的采样点数 30 span = 6; 31 fc = 12e6; %载波频率 32 fs = bit_rate*sps; %采样频率=比特率*每个符号的采样点数 33 rollof_factor = 0.35; %滚降因子,可调整 34 msgLen = 1000; %接收序列长度（包含同步序列） 35 foffset = 30; % Hz 频偏值 36 soffset = 20; % 时偏值 37 deltaderr=0; %时延均方误差 38 deltaferr=0; %频偏均方误差 39 derr=zeros(1,11); %时延均方误差 40 ferr=zeros(1,11); %频偏均方误差 41 42 43 %%%%%%%%%%%%%%%%% 信源 %%%%%%%%%%%%%%%%%%% 44 %%%%本地序列与数据帧的生成 45 %生成本地序列，给出标志性的帧头pn1，方便调试 46 p1=[1,0,0,0,0,0,1]; 47 preamble=mseq(p1);%生成本地序列，周期为127 48 len=length(preamble);%len表示为本地序列的长度，长度定为128 49 if mod(len,2)==1 50 preamble=[preamble 0]; 51 len=len+1; 52 end 53 %将本地序列添加到帧头，生成数据帧 54 playload = randi([0 1],1,872); 55 msg_source= [preamble playload]; 56 57 %%%%%%%%%%%%%%%%% 发射机 %%%%%%%%%%%%%%%%%%% 58 %%%%%调制器 59 %%%双极性变换，就是bpsk映射过程 60 %对I、Q路信号分别进行处理，根据映射关系，始终有Q=0 61 for i=1:msgLen 62 bipolar_msg_source_I(i) = 2*msg_source(i)-1; %双极性信号 63 bipolar_msg_source_Q(i) = 0; 64 end 65 66 %%%%%插值,12倍上采样，原理是在相邻的两个序列的比特之间内插一定数值的0 67 %%%对码元序列进行转置之后在每一行补上相应个数的0 68 bipolar_msg_source_temp_I = [bipolar_msg_source_I',zeros(size(bipolar_msg_source_I,2),sps-1)]; 69 bipolar_msg_source_temp_Q = [bipolar_msg_source_Q',zeros(size(bipolar_msg_source_Q,2),sps-1)]; 70 %%%将补零之后的序列转换一个串行序列 71 length_x_I = size(bipolar_msg_source_temp_I,1); 72 length_y_I = size(bipolar_msg_source_temp_I,2); 73 up12_bipolar_msg_source_I = reshape(bipolar_msg_source_temp_I',1,length_x_I * length_y_I);%长度为12000，12倍插值 74 length_x_Q = size(bipolar_msg_source_temp_Q,1); 75 length_y_Q = size(bipolar_msg_source_temp_Q,2); 76 up12_bipolar_msg_source_Q = reshape(bipolar_msg_source_temp_Q',1,length_x_Q * length_y_Q);%长度为12000，12倍插值 77 78 %%%%%滤波器 79 %%%滚降滤波器 80 rcos_fir = rcosdesign(rollof_factor,span,sps,'sqrt');%长度为73 81 %%%采用滚降滤波器进行滤波 82 rcos_msg_source_I = conv(up12_bipolar_msg_source_I,rcos_fir);%卷积长度为12000+73-1=12072 83 rcos_msg_source_Q = conv(up12_bipolar_msg_source_Q,rcos_fir); 84 rcos_msg_source = rcos_msg_source_I + 1j.*rcos_msg_source_Q ; 85 filter_delay1 = (length(rcos_fir)-1)/2; %滚降滤波器的延迟时长，36 86 %%%生成幅度受到BPSK调制的本地同步序列，并且验证其具有很强的自相关性 87 %采用卷积运算来实现基带成型时，假设数据是一个长度为N的序列a，滤波器是长度为M的序列b，卷积出来的结果长度是N+M-1。 88 %我们要截断出和数据长度N相等的一个长度数据，截断前面的一部分、后面一部分 89 rcos_msg_source_addnoise_prem_down = downsample(rcos_msg_source(filter_delay1+1:end-12*(msgLen-len)-filter_delay1),12);%为降低运算量，选择降采样，12倍，运算量为原来的1/4 90 % corr1=pncorr(rcos_msg_source_addnoise_prem_down,rcos_msg_source_addnoise_prem_down,128); 91 % figure(1) 92 % plot(corr1); 93 94 %%%%%%%%%%%%调制器%%%%%%%%%%%%%%%%%%%% 95 %%%载波发送 96 time = 1:length(rcos_msg_source_I); 97 rcos_msg_source_carrier = rcos_msg_source.*exp(1j*2*pi*fc.*time/fs); 98 99 %%%%%%%%%%%%信道%%%%%%%%%%%%%%%%%%%%%% 100 %%%设置信噪比 101 %SNR(dB)=EbN0(dB)+10log10(K)−10log10(Fs/Rs);%k为几进制，Fs/Rs为过采样倍数 102 ebn0 = 34; 103 snr = ebn0 - 10*log10(sps);%信噪比23dB 104 %%%线性高斯白噪声信道,加入时偏和频偏 105 delay_pad = zeros(1,soffset*sps)+1j*zeros(1,soffset*sps);%手动加入的delay_pad为12*20=240个采样点 106 rcos_msg_source_carrier_sfo_plu = [delay_pad,rcos_msg_source(1:end-soffset*sps)]; 107 rcos_msg_source_carrier_addnoise_plu = awgn(rcos_msg_source_carrier_sfo_plu,snr,'measured');%高斯白噪声信道 108 109 %%%%%%%%%%%%%%接收机%%%%%%%%%%%%%%%%%%%% 110 %%%相干解调,去载波 111 rcos_msg_source_addnoise_plu =rcos_msg_source_carrier_addnoise_plu.*exp(-1j*2*pi*fc.*time/fs); 112 %%%降采样 113 rcos_msg_source_addnoise_down = downsample(rcos_msg_source_addnoise_plu(filter_delay1+1:end-filter_delay1),12); 114 115 %============================================================================ 116 % 生成DPSS参数 117 %%============================================================================ 118 FFO = [-0.05,0.05]; 119 FreqMin = FFO(1); 120 FreqMax = FFO(end); 121 ParW = (FreqMax-FreqMin)/2; 122 ParE = (FreqMax+FreqMin)/2; 123 UseSeq = fix(2*len*ParW)+1; %有效特征值个数2NW+修正；前2NW接近于1，中间有些数，后面基本算0 这里做了改动，加了1 124 [DpssSeq,Lambda] = dpss(len,len*ParW,UseSeq); %问题在这个特征值上面了， 125 %调制DPSS矩阵 126 ColuE = exp(1i*2*pi*ParE*(0:len-1)); %调制矩阵E的对角线元素 127 Dpss = DpssSeq.*repmat(ColuE.',1,UseSeq); 128 EigenValue = real(Lambda); %对角线元素，特征值 129 130 %=================================================== 131 % 时偏估计部分 132 %=================================================== 133 leng_xcorr = msgLen - len; % 相关值的长度（leng_xcorr）= 所有数据的长度 - 同步序列的长度 134 XcorrRes = zeros(1,leng_xcorr); 135 for n= 1:leng_xcorr %滑动相乘 136 for idx = 1:length(EigenValue) %选取的特征值个数 137 LocalSeq = conj(rcos_msg_source_addnoise_prem_down.*Dpss(:,idx).');%乘以权值后的同步序列共轭逆序 138 %corr2=pncorr(LocalSeq,LocalSeq,128); 139 %figure(4) 140 %plot(corr2); 141 %加权之后的同步序列与接收序列做互相关，12个加权之后的同步序列，再将这12个结果相加。 142 xcorr_XX(idx)= abs( rcos_msg_source_addnoise_down((1:len)+n-1)* LocalSeq'); 143 XcorrRes(n)= sum(xcorr_XX);%多特征值合并,累加求和 144 end 145 end 146 %按照分析，应该最大值会出现在20左右， 147 figure(2) 148 plot(XcorrRes);%绘出的是信噪比下的每次移动一个点做互相关运算的结果 149 [max,pos]=max(XcorrRes); 150 151 % %%%直接做互相关运算 152 % for n1= 1:leng_xcorr %滑动相乘 153 % xcorr_X(n1)= abs( rcos_msg_source_addnoise_down((1:len)+n1-1)* conj(rcos_msg_source_addnoise_prem_down).'); 154 % end 155 % figure(5) 156 % plot(xcorr_X) 157 % % 根据仿真结果，可以看出峰值是出现在21，但因为是从1开始取得，这里给的时偏是20，合理的。 158
ok了。时偏估计正确。