二期

  1 clear;
  2 clc;
  3 warning off;
  4 close all;
  5 
  6 length_type=2;
  7 codeandmodtype=10;
  8 data_block=2;%%1-10
  9 Data_punched_length=2;
 10 switch length_type
 11     case 1 ; data_length=1024 ; cp1=128 ;
 12     case 2 ; data_length=1024 ; cp1=256 ;
 13     case 3 ; data_length=1024 ; cp1=512 ;
 14     case 4 ; data_length=2048 ; cp1=128 ;
 15     case 5 ; data_length=2048 ; cp1=256 ;
 16     case 6 ; data_length=2048 ; cp1=512 ;
 17     case 7 ; data_length=4096 ; cp1=128 ;
 18     case 8 ; data_length=4096 ; cp1=256 ;
 19     case 9 ; data_length=4096 ; cp1=512 ;
 20 end
 21 
 22 switch codeandmodtype
 23     case 1 ; code_rate=1/3 ; modtype=1 ;%modtype=1  BPSK
 24     case 2 ; code_rate=1/2 ; modtype=1 ;
 25     case 3 ; code_rate=2/3 ; modtype=1 ;
 26     case 4 ; code_rate=3/4 ; modtype=1 ;
 27     case 5 ; code_rate=4/5 ; modtype=1 ;
 28     case 6 ; code_rate=1/3 ; modtype=2 ;%modtype=2  QPSK
 29     case 7 ; code_rate=1/2 ; modtype=2 ;
 30     case 8 ; code_rate=2/3 ; modtype=2 ;
 31     case 9 ; code_rate=3/4 ; modtype=2 ;
 32     case 10 ; code_rate=4/5 ; modtype=2 ;
 33     case 11 ; code_rate=1/3 ; modtype=3 ;%modtype=3  8PSK
 34     case 12 ; code_rate=1/2 ; modtype=3 ;
 35     case 13 ; code_rate=2/3 ; modtype=3 ;
 36     case 14 ; code_rate=3/4 ; modtype=3 ;
 37     case 15 ; code_rate=4/5 ; modtype=3 ;
 38 end
 39 switch Data_punched_length      %打孔后长度/编码长度
 40     case 1 ; Data_punched_length=2048;   %BPSK / QPSK 所有数据长度
 41     case 2 ; Data_punched_length=4096;   %BPSK 2048及4096 / QPSK所有数据长度
 42     case 3 ; Data_punched_length=6144;   %8PSK 所有数据长度
 43     case 4 ; Data_punched_length=8192;   %BPSK所有数据长度 / QPSK 2048及4096
 44     case 5 ; Data_punched_length=12288;  %8PSK 2048及4096
 45     case 6 ; Data_punched_length=16384;  %QPSK 4096
 46 end
 47 %%  采样速率
 48 in_n = 8;
 49 Fd = 2e6;
 50 Fs = Fd*in_n;   %经过in_n倍过采样,接收到的数据速率是   % Sampling frequency
 51 T = 1/Fs;                 % Sampling period
 52 
 53 %% 对应参数
 54 if cp1==128
 55     zc_length=256 ; cp2=64;
 56 elseif cp1==256
 57     zc_length=512 ; cp2=128;
 58 elseif cp1==512
 59     zc_length=1024 ; cp2=256;
 60 end
 61 
 62 if modtype ==1
 63     hMod =comm.PSKModulator(2, 'BitInput',true,'PhaseOffset',pi/4);%BPSK;
 64 elseif modtype==2
 65     hMod =comm.PSKModulator(4, 'BitInput',true,'PhaseOffset',pi/4);%QPSK;
 66 else
 67     hMod =comm.PSKModulator(8, 'BitInput',true,'PhaseOffset',pi/8);%8PSK;
 68 end
 69 
 70 
 71 %% 莱斯信道
 72 Fsym = Fs;  %% MHz
 73 Ntimes = 1;
 74 fs = Fsym*Ntimes;
 75 RayleighChflag = 1;
 76 offset_flag=1;
 77 AWGNflag = 1;
 78 hRayleighChan_a = comm.RayleighChannel(...
 79     'SampleRate',Fs, ...
 80     'PathDelays',0, ...
 81     'AveragePathGains',0, ...
 82     'NormalizePathGains',true, ...
 83     'MaximumDopplerShift',0, ...
 84     'DopplerSpectrum',doppler('Flat'), ...
 85     'RandomStream','mt19937ar with seed', ...
 86     'Seed',73, ...
 87     'PathGainsOutputPort',false);
 88 hRayleighChan_b = comm.RayleighChannel(...
 89     'SampleRate',Fs, ...
 90     'PathDelays',[0.0 5]*1e-6, ...
 91     'AveragePathGains',[0.0 -6], ...
 92     'NormalizePathGains',true, ...
 93     'MaximumDopplerShift',10, ...
 94     'DopplerSpectrum',doppler('Flat'), ...
 95     'RandomStream','mt19937ar with seed', ...
 96     'Seed',73, ...
 97     'PathGainsOutputPort',false);
 98 hRayleighChan_c = comm.RayleighChannel(...
 99     'SampleRate',Fs, ...
100     'PathDelays',[0 1.5 5.5]*1e-6, ...
101     'AveragePathGains',[0 -5 -9], ...
102     'NormalizePathGains',true, ...
103     'MaximumDopplerShift',25, ...
104     'DopplerSpectrum',doppler('Flat'), ...
105     'RandomStream','mt19937ar with seed', ...
106     'Seed',73, ...
107     'PathGainsOutputPort',false);
108 hRayleighChan_d = comm.RayleighChannel(...
109     'SampleRate',Fs, ...
110     'PathDelays',[0 1.0 2.1 12]*1e-6, ...
111     'AveragePathGains',[0 -2 -5 -9 ], ...
112     'NormalizePathGains',true, ...
113     'MaximumDopplerShift',10, ...
114     'DopplerSpectrum',doppler('Flat'), ...
115     'RandomStream','mt19937ar with seed', ...
116     'Seed',73, ...
117     'PathGainsOutputPort',false);
118 
119 %%  Filter
120 hFilt = comm.RaisedCosineTransmitFilter;   %发送滤波器
121 hFilt.RolloffFactor = 0.4;
122 hFilt.FilterSpanInSymbols = 2;
123 hFilt.OutputSamplesPerSymbol = 8;
124 txCoef=[-0.0007636033507237,-0.001873884698366,-0.002493094497892,-0.002276595999141,-0.001066145936237, 0.001020270601718, 0.003559060766949, 0.005871237463655, 0.007153496030136, 0.006670761767009, 0.003972271435046,-0.0009182805249739, -0.007388879378931, -0.01424081572107, -0.01982327061023, -0.02229012116698,-0.01993954298877, -0.01157441726111, 0.003189250151734,  0.02373440070706,0.04842452872387,  0.07477258410153,  0.09977801998071,0.120378000158,0.133933597037,0.1386619772368, 0.133933597037,  0.120378000158,0.09977801998071,  0.07477258410153,  0.04842452872387,  0.02373440070706,0.003189250151734, -0.01157441726111, -0.01993954298877, -0.02229012116698,-0.01982327061023, -0.01424081572107,-0.007388879378931,-0.0009182805249739,0.003972271435046, 0.006670761767009, 0.007153496030136, 0.005871237463655,0.003559060766949, 0.001020270601718,-0.001066145936237,-0.002276595999141,-0.002493094497892,-0.001873884698366,-0.0007636033507237];
125 
126 RFilt = comm.RaisedCosineReceiveFilter;    %接收滤波器
127 RFilt.RolloffFactor = 0.4;
128 RFilt.FilterSpanInSymbols = 2;
129 RFilt.InputSamplesPerSymbol = 8;
130 RFilt.DecimationFactor = 1;
131 rxCoef=[-0.0007636033507237,-0.001873884698366,-0.002493094497892,-0.002276595999141,-0.001066145936237, 0.001020270601718, 0.003559060766949, 0.005871237463655, 0.007153496030136, 0.006670761767009, 0.003972271435046,-0.0009182805249739, -0.007388879378931, -0.01424081572107, -0.01982327061023, -0.02229012116698,-0.01993954298877, -0.01157441726111, 0.003189250151734,  0.02373440070706,0.04842452872387,  0.07477258410153,  0.09977801998071,0.120378000158,0.133933597037,0.1386619772368, 0.133933597037,  0.120378000158,0.09977801998071,  0.07477258410153,  0.04842452872387,  0.02373440070706,0.003189250151734, -0.01157441726111, -0.01993954298877, -0.02229012116698,-0.01982327061023, -0.01424081572107,-0.007388879378931,-0.0009182805249739,0.003972271435046, 0.006670761767009, 0.007153496030136, 0.005871237463655,0.003559060766949, 0.001020270601718,-0.001066145936237,-0.002276595999141,-0.002493094497892,-0.001873884698366,-0.0007636033507237];
132 
133 
134 % Generate AGC
135 AGC_PN_len = 48;
136 AGC_PN2 =randi([0,1],AGC_PN_len,1);
137 AGC_PN2_mod = pskmod(AGC_PN2, 2, pi/4) * sqrt(2);            %AGC_PN序列BPSK调制  0--> 1 + 1j , 1--> -1 - 1j
138 
139 %% 参数
140 N_sp_code = 16;
141 N_syn_code = 256/N_sp_code;
142 %% 发端
143 r = log2(N_sp_code);
144 n = 1; %第几个多项式
145 poly_nomal=de2bi(primpoly(r,'all'),r+1,'left-msb');
146 ini_state=ones(1,r);
147 pnseq=comm.PNSequence('Polynomial', poly_nomal(n,:), 'SamplesPerFrame', 2^r-1, 'InitialConditions', ini_state);
148 pn_out = [step(pnseq)',0];
149 pn_out_ = pskmod(pn_out, 2, pi/4, 'gray');
150 
151 M1=1;
152 Nzc=16;                                                 %ZC序列长度
153 k=1:Nzc;
154 ZC = zeros(Nzc,1);
155 ZC1(k)=exp((1j*pi*M1*k.^2)/Nzc);
156 sp_code1 = ZC1.';
157 
158 M1=1;
159 Nzc=16;                                                 %ZC序列长度
160 k=[1:Nzc  ];
161 ZC = zeros(Nzc,1);
162 ZC2(1:Nzc)=exp((1j*pi*M1*k.^2)/Nzc);
163 %扩频码
164 sp_code2 = ZC2.';
165 
166 
167 r = log2(N_syn_code);
168 n = 2; %第几个多项式
169 poly_nomal=de2bi(primpoly(r,'all'),r+1,'left-msb');
170 ini_state=ones(1,r);
171 pnseq=comm.PNSequence('Polynomial', poly_nomal(n,:), 'SamplesPerFrame', 2^r-1, 'InitialConditions', ini_state);
172 pn_out = [step(pnseq)',0];
173 data=pn_out;
174 coded_data_m = repmat(data ,length(data),1);
175 coded_data_m = tril(coded_data_m,0);
176 coded_data_x = mod(sum(coded_data_m'),2);
177 pn_out1=1-coded_data_x;
178 pn_out2 = pn_out1*2 - 1;
179 syn_code = pn_out2.';
180 % 同步头
181 syn_head1 = repmat(sp_code1,1,N_syn_code);
182 syn_head1 = syn_head1 .* repmat(syn_code.',N_sp_code,1);
183 syn_head1 = reshape(syn_head1,[],1);
184 
185 
186 syn_head2 = repmat(sp_code2,1,N_syn_code);
187 syn_head2 = syn_head2 .* repmat(syn_code.',N_sp_code,1);
188 syn_head2 = reshape(syn_head2,[],1);
189 
190 % Generate Zadoff-Chu
191 M1=1;
192 Nzc=zc_length;                                                 %ZC序列长度
193 k=1:Nzc;
194 ZC = zeros(Nzc,1);
195 ZC(k)=exp((1j*pi*M1*k.^2)/Nzc);
196 
197 Pilot_PN1=ZC;   %导频序列1
198 Pilot_PNTD1=[Pilot_PN1(end-zc_length/2+1:end);Pilot_PN1];
199 
200 Pilot_PN2=[ZC(end-zc_length/2+1:end);ZC(1:zc_length/2)];   %导频序列1
201 Pilot_PNTD2=[Pilot_PN2(end-zc_length/2+1:end);Pilot_PN2];
202 
203 % Generate head
204 % M1=1;
205 % Nzc=head_length; %     head_length                %ZC序列长度
206 % k=1:Nzc;
207 % ZC_head = zeros(Nzc,1);
208 % ZC_head(k)=exp((1j*pi*M1*k.^2)/Nzc);
209 ZC_head1=Pilot_PNTD1(1:128);
210 ZC_head2=Pilot_PNTD2(1:128);
211 
212 
213 %EbN0= (5:25);
214 
215 % EbN0=[ -2:2:4 5:1:13 13.5:0.5:16];
216 % EbN0= [-4:2:6 7:1:16 16.5:0.5:18];
217 EbN0=[0:1:2 2.5:0.5:3.5 3.7:0.5:10];
218 EbN0=[0:3 3.5:0.5:6];
219 EbN0=[0:1:8 8.5:0.5:12];
220 EbN0=10;
221 %各信噪比仿真次数
222 
223 error_num_snr = zeros(1,length(EbN0));
224 error_num_snr1 = zeros(1,length(EbN0));
225 mu=0;
226 % offset_error = [];
227 for snr_pp =1:length(EbN0)
228     snr_pp;
229     error_num = 0;
230     error_num1 = 0;
231     SNR=zeros(1,length(EbN0));
232     if modtype==1
233         SNR (snr_pp)=EbN0(snr_pp)-10*log10(40/5)-10*log10(1/code_rate);
234     elseif modtype==2
235         SNR (snr_pp)=EbN0(snr_pp)+10*log10(2)-10*log10(40/5)-10*log10(1/code_rate);
236     else
237         SNR (snr_pp)=EbN0(snr_pp)+10*log10(3)-10*log10(40/5)-10*log10(1/code_rate);
238     end
239     SNR;
240     if snr_pp<3
241         n_loop = 10000;
242     elseif snr_pp<9
243         n_loop = 20000;
244     else
245         n_loop = 20000;
246     end
247     n_loop=1000;
248        capture_n1=zeros(1,n_loop);
249     offset_bias = zeros(1,n_loop);
250     f_bias=zeros(1,n_loop);
251     F_begin_max1=zeros(1,n_loop);
252     F_begin_max2=zeros(1,n_loop);
253     F_begin_max11=zeros(1,n_loop);
254     F_begin_max22=zeros(1,n_loop);
255     max_v1=zeros(1,n_loop);
256     max_v2=zeros(1,n_loop);
257     
258     eee=0;
259     for k = 1:n_loop
260     
261         [Data_punched1,codeword_order,data_org]=codemodtype(data_length,data_block,Data_punched_length,modtype,code_rate);
262         data_moded =step(hMod, reshape(Data_punched1,[],1));   % mod
263         data_mod=reshape(data_moded,[],data_block);
264         
265         
266         %% 多天线
267         d1=data_mod(1:data_length,:);
268         d2=data_mod(data_length+1:end,:);
269         
270         d1_conj=conj(d1);
271         d2_conj=-conj(d2);
272         d1_conj_daoxu=flip(d1_conj);
273         d1_conj_daoxu=[d1_conj_daoxu(end,:);d1_conj_daoxu(1:end-1,:)];
274         d2_conj_daoxu=flip(d2_conj);
275         d2_conj_daoxu=[d2_conj_daoxu(end,:);d2_conj_daoxu(1:end-1,:)];
276         tx_1=[d1(end-cp2+1:end,:);d1;d2_conj_daoxu(end-cp2+1:end,:);d2_conj_daoxu];
277         tx_2=[d2(end-cp2+1:end,:);d2;d1_conj_daoxu(end-cp2+1:end,:);d1_conj_daoxu];
278         %         SentSignal1 = [AGC_PN2_mod;Pilot_PNTD1;data_mod(1:2176,:);Pilot_PNTD1;data_mod(2177:end,:)];
279         %% construct frame Tx1
280         if data_block==1
281             data_tx1=[AGC_PN2_mod;(syn_head1);Pilot_PNTD1;tx_1;Pilot_PNTD1];
282         else
283             pilot_add_data1=zeros(data_length*2+cp2*2+zc_length*1.5,data_block);
284             for aa=1:data_block
285                 pilot_add_data1(:,aa)=[tx_1(:,aa);Pilot_PNTD1];
286             end
287             pilot_add_data1=reshape(pilot_add_data1,[],1);
288             data_tx1=[AGC_PN2_mod;(syn_head1);Pilot_PNTD1;pilot_add_data1];
289         end
290         up_data_tx1=upsample(data_tx1,8,0);
291         SentSignal_1=conv(txCoef,up_data_tx1);
292         SentSignal1 = SentSignal_1(26:end-25);
293         %         SentSignal = step(hFilt,SentSignal1);
294         
295         
296         %% construct frame Tx2
297         if data_block==1
298             data_tx2=[AGC_PN2_mod;(syn_head2);Pilot_PNTD2;tx_2;Pilot_PNTD2];
299         else
300             pilot_add_data2=zeros(data_length*2+cp2*2+zc_length*1.5,data_block);
301             for aa=1:data_block
302                 pilot_add_data2(:,aa)=[tx_2(:,aa);Pilot_PNTD2];
303             end
304             pilot_add_data2=reshape(pilot_add_data2,[],1);
305             data_tx2=[AGC_PN2_mod;(syn_head2);Pilot_PNTD2;pilot_add_data2];
306         end
307         up_data_tx2=upsample(data_tx2,8,0);
308         SentSignal_2=conv(txCoef,up_data_tx2);
309         SentSignal2 = SentSignal_2(26:end-25);
310         %% Create an AWGN channel.
311         %%%% 加频偏
312         offset_tx1 =8000;
313         offset_tx2 =8000;
314         if offset_flag == 1
315             
316             Signal_through_channel_add_offset1 = SentSignal1.' .* exp( (1:1:length(SentSignal1)) * 1j * 2 * pi * offset_tx1 / Fs);
317             %         Signal_through_channel_add_offset1=[zeros(1,1000), Signal_through_channel_add_offset1];
318             
319             Signal_through_channel_add_offset2 = SentSignal2 .'.* exp( (1:1:length(SentSignal2)) * 1j * 2 * pi * offset_tx2 / Fs);
320             %         Signal_through_channel_add_offset2=[zeros(1,1000), Signal_through_channel_add_offset2];
321         else
322             Signal_through_channel_add_offset1=SentSignal1.';
323             Signal_through_channel_add_offset2=SentSignal2.';
324         end
325         %%%  加莱斯信道
326         if RayleighChflag == 1
327             Sig_tx_01 = step(hRayleighChan_b,Signal_through_channel_add_offset1.');
328             Sig_tx_02 = step(hRayleighChan_b,Signal_through_channel_add_offset2.');
329             Sig_tx_01=Sig_tx_01.';
330             Sig_tx_02=Sig_tx_02.';
331         else
332             Sig_tx_01 = Signal_through_channel_add_offset1;
333             Sig_tx_02 = Signal_through_channel_add_offset2;
334         end
335         
336         %%%%  加高斯噪声
337         %         sigpower1=mean(abs(Sig_tx_01(48*8+128*8+1:end)).^2);
338         %         sigpower1=10*log10(sigpower1);
339         if AWGNflag == 1
340             Signal_through_channel_noised_add_offset1=awgn(Sig_tx_01,SNR(snr_pp),'measured');
341             Signal_through_channel_noised_add_offset2=awgn(Sig_tx_02,SNR(snr_pp),'measured');
342             %              Signal_through_channel_noised_add_offset=Sig_tx_01.';
343         else
344             Signal_through_channel_noised_add_offset1=Sig_tx_01;
345             Signal_through_channel_noised_add_offset2=Sig_tx_02;
346         end
347         
348         %% 匹配滤波
349         Signal_after_matched_filter_1=conv( rxCoef,Signal_through_channel_noised_add_offset1);
350         Signal_after_matched_filter1 = Signal_after_matched_filter_1(26:end-25);
351         %         Signal_after_AD1 = Signal_after_matched_filter1.';
352         
353         Signal_after_matched_filter_2=conv( rxCoef,Signal_through_channel_noised_add_offset2);
354         Signal_after_matched_filter2 = Signal_after_matched_filter_2(26:end-25);
355         %         Signal_after_AD2 = Signal_after_matched_filter2.';
356         
357         Signal_after_matched_filter=Signal_after_matched_filter1+Signal_after_matched_filter2;
358         %         Signal_after_matched_filter=Signal_after_matched_filter1;
359         %% 抽取为4倍
360         Signal_after_matched_filter_4up=Signal_after_matched_filter(1:2:end);
361         Signal_after_matched_filter_4up=Signal_after_matched_filter_4up.';
362         in_n1=4;
363        %% 8相位量化  不仅影响同步性能 而且影响ber  损失至少1dB
364         rx = Signal_after_matched_filter_4up;
365         angle_rx = angle(rx);
366         
367         rx(find(angle_rx<pi/8 & angle_rx>=-pi/8)) = 1;
368         rx(find(angle_rx<3*pi/8 & angle_rx>=pi/8)) = (1 + 1i)/sqrt(2);
369         rx(find(angle_rx<5*pi/8 & angle_rx>=3*pi/8)) =  1i;
370         rx(find(angle_rx<7*pi/8 & angle_rx>=5*pi/8)) =  (-1 + 1i)/sqrt(2);
371         rx(find(-angle_rx>7*pi/8 | angle_rx>=7*pi/8)) =  -1;
372         rx(find(-angle_rx<=3*pi/8 & -angle_rx>pi/8)) = (1 - 1i)/sqrt(2);
373         rx(find(-angle_rx<=5*pi/8 & -angle_rx>3*pi/8)) =  -1i;
374         rx(find(-angle_rx<=7*pi/8 & -angle_rx>5*pi/8)) =  (-1 - 1i)/sqrt(2);
375         Signal_rx_filt_LH = rx;
376         
377         
378         %% ZC1差分相关
379         %帧同步
380         %解扩 第一次相关
381         First_corr1 =[];
382         for mm = 1:length(Signal_rx_filt_LH)-(N_sp_code-1)*in_n1
383             First_corr1(mm) = sum(Signal_rx_filt_LH(mm:in_n1:mm+(N_sp_code-1)*in_n1).*conj(sp_code1));
384             %     First_corr_power(k) = sqrt(sum(abs(Signal_rx_filt_LH(k:in_n:k+(N_sp_code-1)*in_n)).^2));
385         end      
386 
387         First_corr_abs = abs(First_corr1).^2;
388          
389         % figure, stem(First_corr_abs);
390         % figure, plot(First_corr,'*');
391         % 第二次相关
392         
393 %         rx = First_corr1;
394 %         angle_rx = angle(rx);
395 %         % 8相位量化  不仅影响同步性能 而且影响ber  损失至少1dB
396 %         rx(find(angle_rx<pi/8 & angle_rx>=-pi/8)) = 1;
397 %         rx(find(angle_rx<3*pi/8 & angle_rx>=pi/8)) = (1 + 1i)/sqrt(2);
398 %         rx(find(angle_rx<5*pi/8 & angle_rx>=3*pi/8)) =  1i;
399 %         rx(find(angle_rx<7*pi/8 & angle_rx>=5*pi/8)) =  (-1 + 1i)/sqrt(2);
400 %         rx(find(-angle_rx>7*pi/8 | angle_rx>=7*pi/8)) =  -1;
401 %         rx(find(-angle_rx<=3*pi/8 & -angle_rx>pi/8)) = (1 - 1i)/sqrt(2);
402 %         rx(find(-angle_rx<=5*pi/8 & -angle_rx>3*pi/8)) =  -1i;
403 %         rx(find(-angle_rx<=7*pi/8 & -angle_rx>5*pi/8)) =  (-1 - 1i)/sqrt(2);
404 % %         First_corr1 = rx;
405 %             First_corr1 = First_corr1;
406 %         First_corr1_abs = abs(First_corr1).^2;
407         First_corrx1 = First_corr1(1:end-in_n1*N_sp_code).*conj(First_corr1(1+in_n1*N_sp_code:end));
408         
409         Second_xcorr1=[];
410             syn_code_x = 2*pn_out-1;
411         for mm = 1:length(First_corrx1)-(N_syn_code-1)*in_n1*N_sp_code
412             xxx = First_corrx1(mm:in_n1*N_sp_code:mm+(N_syn_code-1)*in_n1*N_sp_code);
413             Second_xcorr1(mm) = abs(sum(xxx.*(syn_code_x))).^2;
414         end
415         % figure, stem(Second_xcorr1);
416 %         Second_xcorr1=[zeros(1,N_sp_code*in_n1),Second_xcorr1];
417         [max_v1(k),F_begin_max1(k)] =max(Second_xcorr1);
418 %         for ii=1:4*cp2
419 %             if Second_xcorr1(F_begin_max11)
420 %           [max_v1(k), F_begin_max1(k)]=  max(Second_xcorr1(F_begin_max11-4*cp2:F_begin_max11)) ;
421         syn_n1 = F_begin_max1(k)+N_sp_code*in_n1;
422 %         syn_n1=syn_n1-cp2*2;
423 %         syn_n1=193;
424         
425         if abs(syn_n1-193)>2 && abs(syn_n1-545)>2
426             eee=eee+1;
427 %             disp(['  syn_n: ', num2str(syn_n1)]);
428         end
429         % 计算频偏
430         corr_peak=First_corr1(syn_n1:N_sp_code*in_n1:N_sp_code*in_n1*(N_syn_code-1)+syn_n1).*syn_code';
431 %         dlmwrite('corr_peak.txt',corr_peak,'delimiter','');
432         fft_N = 2048;
433         [corr_fre_est1, fft_n_fre_est1] = max(abs(fft(First_corr1(syn_n1:N_sp_code*in_n1:N_sp_code*in_n1*(N_syn_code-1)+syn_n1).*syn_code',fft_N)));%.*syn_code'
434 
435         [max_fft,num1] = max(corr_fre_est1);
436         fft_n1=fft_n_fre_est1(num1);
437         
438         if fft_n1>fft_N/2
439             fre_set1 = -(fft_N - fft_n1+1)*Fd/N_sp_code/fft_N;%计算频偏
440         else
441             fre_set1 = -(fft_n1-1)*Fd/N_sp_code/fft_N;
442         end
443         timing_loc1=syn_n1+N_sp_code*in_n1*N_syn_code;
444 %         f_bias(k)=fre_set1+offset_tx1;
445 %% ZC2差分相关
446         %帧同步
447         %解扩 第一次相关
448         First_corr2 =[];
449         for mm = 1:length(Signal_rx_filt_LH)-(N_sp_code-1)*in_n1
450             First_corr2(mm) = sum(Signal_rx_filt_LH(mm:in_n1:mm+(N_sp_code-1)*in_n1).*conj(sp_code2));
451             %     First_corr_power(k) = sqrt(sum(abs(Signal_rx_filt_LH(k:in_n:k+(N_sp_code-1)*in_n)).^2));
452         end      
453         First_corr2=First_corr2;
454         First_corr_abs2 = abs(First_corr2).^2;
455          
456         % figure, stem(First_corr_abs2);
457         % figure, plot(First_corr,'*');
458         % 第二次相关
459         
460         rx = First_corr2;
461         angle_rx = angle(rx);
462         % 8相位量化  不仅影响同步性能 而且影响ber  损失至少1dB
463         rx(find(angle_rx<pi/8 & angle_rx>=-pi/8)) = 1;
464         rx(find(angle_rx<3*pi/8 & angle_rx>=pi/8)) = (1 + 1i)/sqrt(2);
465         rx(find(angle_rx<5*pi/8 & angle_rx>=3*pi/8)) =  1i;
466         rx(find(angle_rx<7*pi/8 & angle_rx>=5*pi/8)) =  (-1 + 1i)/sqrt(2);
467         rx(find(-angle_rx>7*pi/8 | angle_rx>=7*pi/8)) =  -1;
468         rx(find(-angle_rx<=3*pi/8 & -angle_rx>pi/8)) = (1 - 1i)/sqrt(2);
469         rx(find(-angle_rx<=5*pi/8 & -angle_rx>3*pi/8)) =  -1i;
470         rx(find(-angle_rx<=7*pi/8 & -angle_rx>5*pi/8)) =  (-1 - 1i)/sqrt(2);
471 %         First_corr1 = rx;
472             First_corr2 = First_corr2;
473         First_corr2_abs = abs(First_corr2).^2;
474         First_corrx2 = First_corr2(1:end-in_n1*N_sp_code).*conj(First_corr2(1+in_n1*N_sp_code:end));
475         
476         Second_xcorr2=[];
477             syn_code_x = 2*pn_out-1;
478         for mm = 1:length(First_corrx2)-(N_syn_code-1)*in_n1*N_sp_code
479             xxx = First_corrx2(mm:in_n1*N_sp_code:mm+(N_syn_code-1)*in_n1*N_sp_code);
480             Second_xcorr2(mm) = abs(sum(xxx.*(syn_code_x))).^2;
481         end
482         % figure, stem(Second_xcorr2);
483         % figure, stem(power_normal);
484         % figure, stem(Second_xcorr_normal);
485         
486         [max_v2(k),F_begin_max2(k)] =max(Second_xcorr2);
487 %          F_begin_max2(k)=  max(Second_xcorr2(F_begin_max22-4*cp2+N_sp_code*in_n1:F_begin_max22)) ;
488         syn_n2 = F_begin_max2(k)+N_sp_code*in_n1;
489 %         syn_n2=syn_n2-cp2*2;
490 %         syn_n2=193;
491         
492 %         if abs(syn_n2-193)>2 && abs(syn_n2-545)>2
493 %             eee=eee+1;
494 % %             disp(['  syn_n: ', num2str(syn_n2)]);
495 %         end
496         % 计算频偏
497         
498         fft_N = 2048;
499         [corr_fre_est2, fft_n_fre_est2] = max(abs(fft(First_corr2(syn_n2:N_sp_code*in_n1:N_sp_code*in_n1*(N_syn_code-1)+syn_n2).*syn_code',fft_N)));%.*syn_code'
500 
501         [max_fft,num2] = max(corr_fre_est2);
502         fft_n2=fft_n_fre_est2(num2);
503         
504         if fft_n2>fft_N/2
505             fre_set2 = -(fft_N - fft_n2+1)*Fd/N_sp_code/fft_N;%计算频偏
506         else
507             fre_set2 = -(fft_n2-1)*Fd/N_sp_code/fft_N;
508         end
509         timing_loc2=syn_n2+N_sp_code*in_n1*N_syn_code;
510 %         f_bias(k)=fre_set1+offset_tx1;
511 fre_set=(fre_set1+fre_set2)/2;
512 
513               timing_loc=min(timing_loc1,timing_loc2);
514 %               timing_loc=1217;
515 %        Signal_after_matched_filter_4up1=[Signal_after_matched_filter_4up;zeros(1000,1)] ; 
516 %         Signal_after_matched_filter_4up1=Signal_after_matched_filter_4up1(timing_loc:end);
517 %        max_v1(k)
518 %        max_v2(k)
519         %% 下采
520         Signal_after_matched_filter_4up1=[Signal_after_matched_filter_4up;zeros(1000,1)];
521 %         Signal_after_matched_filter_4up=Signal_after_matched_filter_4up1(timing_loc:end);
522         
523 
524      
525             down_Signal_after_AD2=Signal_after_matched_filter_4up1(timing_loc:2:timing_loc+(cp1+zc_length)*(data_block+1)*in_n1+(cp2+data_length)*data_block*2*in_n1-1);
526             %往前提0.5CP
527 %                  down_Signal_after_AD2=Signal_after_matched_filter_4up1(timing_loc-cp2*2:2:timing_loc-cp2*2+(cp1+zc_length)*(data_block+1)*in_n1+(cp2+data_length)*data_block*2*in_n1-1);
528              %定时定准
529 %                  down_Signal_after_AD2=Signal_after_matched_filter_4up1(1217:2:1217+(cp1+zc_length)*(data_block+1)*in_n1+(cp2+data_length)*data_block*2*in_n1-1);
530       
531 
532  %% 粗纠(频域估计的频偏)
533         Signal_after_roughFTR=down_Signal_after_AD2.'.* exp( (1:1:length(down_Signal_after_AD2)) * 1j * 2 * pi * fre_set / Fs*4);           
534                exact_total_Signal=Signal_after_roughFTR;
535        %% 粗频偏估计(分别用block_num+1个导频的cp各自估计的频偏作为夹在中间的数据的频偏)
536         in_n1=2;
537         rough_esti_corr=zeros(1,data_block);
538         rough_offset_esti=zeros(1,data_block);
539         LL2=zc_length*in_n1;
540         rough_correct_len=zc_length*in_n1/2;
541         zc_interval1=(zc_length*1.5+(cp2+data_length)*2)*in_n1;
542         num_2=1;
543         last_pilot_pos=length(exact_total_Signal);
544         for i=1:zc_interval1:last_pilot_pos
545             rough_esti_zc1=exact_total_Signal(i: i+ rough_correct_len - 1);                 %用接收序列相关进行频偏估计取相邻1024
546             rough_esti_zc2=exact_total_Signal(LL2+i:LL2+i+ rough_correct_len-1);
547             rough_esti_corr( num_2)=sum(rough_esti_zc1.*conj(rough_esti_zc2))*exp(-1i*0/Fs *2*LL2*2*pi);
548             num_2= num_2+1;
549         end
550         rough_esti_angle1=angle(rough_esti_corr);
551         rough_offset_esti1 = rough_esti_angle1 * Fs /4/ 2 / pi ./ LL2;
552         
553         for i=1:data_block
554             rough_offset_esti(i)=mean(rough_offset_esti1(i:i+1));
555         end
556 
557         %         offset_error1 = rough_offset_esti+offset;
558         %粗纠
559 %         Signal_after_roughFTR = Signal_after_matched_filter_4up.'.*exp((1:1:length(Signal_after_matched_filter_4up))*1j*2*pi*(0)/Fs*2);
560        %% 精频偏估计(用两段导频估计)
561         exact_esti_corr=zeros(1,data_block);
562         LL3=(cp1+zc_length+(cp2+data_length)*2)*in_n1;
563         exact_correct_len=zc_length*in_n1;
564         zc_interval2=(zc_length*1.5+(cp2+data_length)*2)*in_n1;
565         num_3=1;
566         last_pilot_pos=length(exact_total_Signal) - zc_length*1.5*in_n1;
567         for i=1:data_block
568             exact_esti_zc1=exact_total_Signal(zc_length*in_n1/2+1+zc_interval2*(i-1):zc_length*in_n1/2+exact_correct_len+zc_interval2*(i-1));
569             exact_esti_zc2=exact_total_Signal(zc_length*in_n1/2+1+zc_interval2*(i-1)+LL3:zc_length*in_n1/2+exact_correct_len+LL3+zc_interval2*(i-1));
570             exact_esti_corr( num_3)=sum(exact_esti_zc1.*conj(exact_esti_zc2))*exp(-1i*(rough_offset_esti(i))/Fs *4*LL3*2*pi);
571            num_3= num_3+1;
572         end
573 
574         exact_esti_angle=angle(exact_esti_corr);
575         exact_offset_esti = exact_esti_angle * Fs /4/ 2 / pi ./ LL3;
576                 
577 
578         %% 细纠
579         offset_total_exact=zeros(1,data_block);
580         Signal_after_FTR=zeros(length(down_Signal_after_AD2),1);
581         LL4=(cp1+zc_length+cp2*2+data_length*2)*in_n1;
582         down_Signal_after_AD1=reshape(down_Signal_after_AD2(1:end-(cp1+zc_length)*2),[],data_block);
583         for i=1:data_block
584             offset_total_exact(i)=exact_offset_esti(i) + rough_offset_esti(i) ;
585         end
586         num_3=1;
587         for i=1:LL4:length(down_Signal_after_AD2)-(cp1+zc_length)*2
588             if i==(data_block-1)*LL4+1
589                 signal_temp= exact_total_Signal(i:i+LL4-1+(cp1+zc_length)*2);
590                 Signal_after_FTR(i:i+LL4-1+(cp1+zc_length)*2) = signal_temp .* exp( (i:1:i+length(signal_temp)-1) * 1j * 2 * pi * offset_total_exact(num_3) / Fs*4);            %细纠
591             else 
592                 signal_temp= exact_total_Signal(i:i+LL4-1);
593                 Signal_after_FTR(i:i+LL4-1) = signal_temp .* exp( (i:1:i+length(signal_temp)-1) * 1j * 2 * pi * offset_total_exact(num_3) / Fs*4);            %细纠
594             end
595             num_3=num_3+1;
596         end
597         down_Signal_after_AD_after_FTR=[zeros(cp2,1); Signal_after_FTR(1:end-cp2)];
598         
599             
600        
601         %         offset_error(snr_pp,k) = rough_offset_esti+exact_offset_esti+offset;
602 %         offset_bias(k)=rough_offset_esti+exact_offset_esti+offset;
603    
604         
605         %% 取出数据和导频 去cp
606         pilot_receive_all = zeros((cp1+zc_length)*2,data_block+1);
607         data_receive_all = zeros((cp2+data_length)*2*2,data_block);
608         
609         down_Signal_after_AD_after_FTR=[zeros(cp2,1);down_Signal_after_AD_after_FTR(cp2+1:end)];
610         down_Signal_after_AD_after_FTR_11=reshape(down_Signal_after_AD_after_FTR(1:end-(cp1+zc_length)*2),[],data_block);
611         
612         for ii=1:data_block
613             pilot_receive_all(:,ii)=down_Signal_after_AD_after_FTR_11(1:(cp1+zc_length)*2,ii);
614             data_receive_all(:,ii)=down_Signal_after_AD_after_FTR_11(((cp1+zc_length)*2+1:end),ii);%(cp1+zc_length)*2+(cp1+data_length)*4
615         end
616         pilot_receive_all(:,data_block+1)=down_Signal_after_AD_after_FTR(end-(cp1+zc_length)*2+1:end);
617         
618         pilot_receive=pilot_receive_all(cp1*2+1:end,:);
619         data_receive_cp=reshape(data_receive_all,[],data_block*2);
620         data_receive=data_receive_cp(cp2*2+1:end,:);
621         
622         for i=1:data_block
623             data_receive1=data_receive(:,1:2:end);
624             data_receive2=data_receive(:,2:2:end);
625         end
626         
627         
628         %% Channel estimation 收端两段导频(去cp)和本地zc序列相除得出的值作为夹着的数据的信道估计值(LS估计)
629         
630         pilot_demod=zeros(zc_length*2,data_block+1);
631         pilot_demod_zclength=zeros(zc_length,data_block+1);
632         H_esti1_zclength1=zeros(zc_length,data_block+1);
633         H_esti1_zclength=zeros(zc_length,data_block);
634         FXH1=zeros(zc_length,data_block+1);
635         for i=1:data_block+1
636             pilot_receive_temp =fft(pilot_receive(:,i))./sqrt(2*zc_length);%    %plot(abs(pilot_receive_temp))
637             pilot_demod(:,i)=pilot_receive_temp;
638             %     pilot_demod_256(:,i)=[pilot_demod(129:128+256,i)];%
639             pilot_demod_zclength(:,i)=[pilot_demod(1:zc_length/2,i);pilot_demod(zc_length*1.5+1:end,i)];%==时域滤波降采
640         end
641         
642         %         FXH21=pilot_demod_2561(:,2);
643         ZC256FFT=fft(Pilot_PN1)./sqrt(zc_length);
644         for ii=1:data_block+1
645             FXH1(:,ii)=pilot_demod_zclength(:,ii);
646             H_esti1_zclength1(:,ii)=FXH1(:,ii)./(ZC256FFT);%LS估计
647         end
648         for iii=1:data_block
649             H_esti1_zclength(:,iii)=(H_esti1_zclength1(:,iii)+H_esti1_zclength1(:,iii+1))/2;
650         end
651 %                 H_esti1_zclength=H_esti1_zclength1(:,1:data_block);
652         H_esti_2zclength=[H_esti1_zclength(1:zc_length/2,:);zeros(zc_length,data_block);H_esti1_zclength(zc_length/2+1:end,:)];%==时域上采滤波
653         h_esti=ifft(H_esti_2zclength).*sqrt(zc_length*2);
654         shift_n = zc_length/2-0;
655         h_esti1=[h_esti(1:zc_length-shift_n,:);zeros(data_length*2-zc_length,data_block);h_esti(end-shift_n+1:end,:)];%12,00,78
656         h_esti2=[h_esti(zc_length+1:end-shift_n,:);zeros(data_length*2-zc_length,data_block);h_esti(zc_length+1-shift_n:zc_length,:)];%56,00,34
657         
658         n_var=(sum(abs(h_esti(zc_length/2+1:zc_length*0.75)).^2)+sum(abs(h_esti(zc_length*1.5+1:zc_length*1.75)).^2))/(zc_length/2);%3,7
659         threshold=n_var*8;
660         h_esti1(abs(h_esti1).^2<threshold)=0;
661         h_esti2(abs(h_esti2).^2<threshold)=0;
662         
663 %          h_esti_power1=(mean(abs(h_esti1).^2)+mean(abs(h_esti2).^2));
664 %                         h_esti_power=(mean(abs(h_esti).^2));
665 %                         data_receive_pow=(mean(abs(data_receive).^2));
666 %                 SNR_esti1=h_esti_power1/n_var;
667 %                 SNR_esti=10*log10(SNR_esti1);
668 %                 SNR;  %信噪比估计
669 %          dlmwrite('h_esti1_i.txt',real(h_esti1(:,1)),'delimiter','');
670 %          dlmwrite('h_esti1_q.txt',imag(h_esti1(:,1)),'delimiter','');
671 %          h_esti1_i=load('H:\研\单载波\双发\4.12加编解码+ZC16扩频码\h_esti1_i.txt');
672 %          h_esti1_q=load('H:\研\单载波\双发\4.12加编解码+ZC16扩频码\h_esti1_q.txt');
673 %          h_esti1(:,1)=h_esti1_i+1i*h_esti1_q;
674         H_esti11_1=fft(h_esti1)/sqrt(length(h_esti1));
675         H_esti11_2=fft(h_esti2)/sqrt(length(h_esti1));
676         
677         %         H_esti11_pow=mean(abs(H_esti11_1).^2);
678         %         h_esti_pow=mean(abs(h_esti).^2);
679         %         aaaa=H_esti11_pow/h_esti_pow;
680         
681         %         H_esti11_11=repmat(H_esti11_1,1,2);
682         
683         %% Y1、Y2处理
684         data_demod1=zeros(data_length*2,data_block);
685         data_demod2=zeros(data_length*2,data_block);
686         for i=1:data_block
687             data_receive_temp1 =fft(data_receive1(:,i))./sqrt(data_length*2);
688             %     data_receive_temp=fftshift(data_receive_temp);
689             data_demod1(:,i)=data_receive_temp1;
690         end
691         
692         for i=1:data_block
693             data_receive2_daoxu=flip(data_receive2(:,i));
694             %             data_receive2_daoxu=[data_receive2_daoxu(end,:);data_receive2_daoxu(1:end-1,:)];
695             data_receive2_daoxu_c=conj(data_receive2_daoxu);
696             data_receive2_daoxu_c=[data_receive2_daoxu_c(end,:);data_receive2_daoxu_c(1:end-1,:)];
697             data_receive_temp2 =fft(data_receive2_daoxu_c)./sqrt(data_length*2);
698             %             data_receive_temp2=[data_receive_temp2(end,:);data_receive_temp2(1:end-1,:)];
699             %     data_receive_temp=fftshift(data_receive_temp);
700             data_demod2(:,i)=data_receive_temp2;
701         end
702         
703         
704         %% 分离D1与D2
705         H_esti11_1_c=conj(H_esti11_1);
706         H_esti11_2_c=conj(H_esti11_2);
707         
708         R1=data_demod1.*H_esti11_1_c + data_demod2.*H_esti11_2;%均衡
709         R2=data_demod1.*H_esti11_2_c - data_demod2.*H_esti11_1;
710         
711         Habs_1=H_esti11_1.*H_esti11_1_c+n_var/4;
712         Habs_2=H_esti11_2.*H_esti11_2_c+n_var/4;
713         D=zeros(data_length*2,data_block*2);
714         D(:,1:2:end)=R1./(Habs_1+Habs_2);
715         D(:,2:2:end)=R2./(Habs_1+Habs_2);
716         D_1024=[D(1:data_length/2,:);D(data_length*1.5+1:end,:)];
717         %          D_r=[D_1024(:,1),D_1024(:,2)];
718         D_shiyu=ifft(D_1024).*sqrt(data_length);
719         
720         %          R1=data_demod1 + data_demod2;%
721         %          R2=data_demod1 - data_demod2;
722         %          D2222=zeros(2048,2);
723         %          D2222(:,1)=R1./2;
724         %          D2222(:,2)=R2./2;
725         %          D_102422=[D2222(1:512,:);D2222(1537:end,:)];
726         % %          D_r=[D_1024(:,1),D_1024(:,2)];
727         %          D_shiyu22=ifft(D_102422).*sqrt(102422);
728         
729         
730         %% QPSK modulation
731         %         %%%%
732         %         data_receive1=data_receive(1:2:end);
733         %                 data_de_ploit=reshape(data_receive1,[],1);   % 无均衡
734         %                     scatterplot(reshape(D_shiyu,[],1));
735         data_de_ploit=reshape(D_shiyu,[],1); %有均衡
736         %         data_decoded=zeros(2048,3);
737         
738         
739         
740         %软判
741         if modtype==1
742             %             data_demod_bit=zeros(data_length*2*data_block,1);
743             data_demod_bit=real(data_de_ploit)+imag(data_de_ploit);
744         elseif modtype==2
745             data_demod_bit=zeros(data_length*4*data_block,1);
746             for llr_i=1:length(data_de_ploit)
747                 data_demod_bit(llr_i*2,:)=real(data_de_ploit(llr_i));
748                 data_demod_bit(llr_i*2-1,:)=imag(data_de_ploit(llr_i));
749             end
750             data_demod_bit=reshape(data_demod_bit,[],data_block);
751         elseif modtype==3
752             niose_var=n_var/4;
753             hDemod = comm.PSKDemodulator(8,'PhaseOffset',pi/8, 'BitOutput',true,...
754                 'DecisionMethod','Approximate log-likelihood ratio', ...
755                 'Variance', niose_var);%8PSK;
756             data_demod_bit = step(hDemod,data_de_ploit);
757         end
758         
759         %&解码
760 
761             [rx_bit]=decode(data_length,data_block,modtype,code_rate,Data_punched_length,data_demod_bit);
762        total_data=length(reshape(rx_bit,[],1));
763         rx_bit=reshape(rx_bit,[],1);
764         error_num(k) = sum(abs(rx_bit~=data_org));
765 %                 dbstop in danzaibo_2tsyntiqian at 776 if error_num(k)~=0;
766         
767         %     toc;
768     end
769 %     max_v1
770 %        max_v2
771     %     delete(gcp('nocreate'))
772     capture_nmax=max(capture_n1);
773     capture_nmin=min(capture_n1); 
774     error_num_snr(snr_pp) = sum(error_num);
775     error_rate(snr_pp)=error_num_snr(snr_pp)/n_loop/total_data;
776 end
777 
778 %  data_length=13824;
779 
780 semilogy(EbN0 , error_rate,'r-^');hold on
781 xlabel('EbN0');
782 ylabel('error');
783 title('data=1024 zc=256 QPSK 0.8_4096精纠和均衡定时定准');
784 grid on;

 

上一篇:实验2 数组、指针与c++标准库


下一篇:shell中的单括号[ ]、 双括号[[ ]] 和 test的区别