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;