北斗信号捕获

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
%载波 rom
Ac = 1024;
c = 0:1/1023:1;
c_rom = round(Ac*cos(2*pi*c));
s_rom = round(Ac*cos(2*pi*c));
CwFreq = CodeFreq;%载波频率；
NHCode = [0 0 0 0 0 1 0 0 1 1 0 1 0 1 0 0 1 1 1 0];%NH码
NHCode(NHCode == 0) = -1;%NH码 为0时，就变为-1，
NHCode = repelem(NHCode,1,DataLen); %按照数据长度进行数据扩展，为什么要扩展呢，原因是我们生成的BDData  1行DataLen列
%要与我们的NH码进行相乘，而BDData数据长度是10，所以要将NH码进行数据数据长度的扩展，
BDData = rand(1,DataLen);%随机产生10个帧数据，，需要转换到-1,1；
BDData(BDData < 0.5) = -1;%将小于0.5的数据都变为-1，
BDData(BDData > 0.5) = 1;%将大于0.5的数据都变为1，
BDData = repelem(BDData, NHLen);
%NH码调制过程
ModNH = NHCode.*BDData; %NHLen*DataLen
ModNH = repelem(ModNH, (SampleFreq/CodeFreq)*CodeLength); %升采样，与后Pn1长度相等
%产生扩频码，为什么是用那个函数产生
SignalLen = (SampleFreq/CodeFreq)*CodeLength*NHLen*DataLen;%总的数据长度，信号的中长度满找采样率产生
[pn1,pn1_freq] = pn_code_gen(PRN_num1,SignalLen,SampleFreq,CodeFreq,CodeLength); %PRN6
ModPN = ModNH.*pn1_freq; %把NH调制之后的数据进行PN调制
% S =  A*C*D*cos(2*pi*fc + )  C:PN码，伪随机码，乘以NH,D：导航数据，进行BPSK,调制
%频率控制字
CWFreq_word  = CwFreq * 2^32 /SampleFreq; %频率控制字
dac_o = zeros(1,SignalLen);
wr = 0;
romaddr = 1;
for i = 1:SignalLen
    if(ModPN(i) == 1)
        dac_o(i) = c_rom(romaddr); %BPSk调制
    elseif(ModPN(i) == -1)
        dac_o(i) = -c_rom(romaddr);
    end
    wr = wr + CWFreq_word;
    if(wr > 2^32)
        wr = wr - 2^32;
    end
    romaddr = round((wr/2^22));
    if(romaddr == 0)
        romaddr = 1; 
    end
end
%%%%%%%%读文件
file = fopen('TEST6.BIN','rb');
data = fread(file,'int16');
fclose(file);
dac_data = data';%%转置，1行N列，转为1列N行

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
%数据偏移
dataoffset = 2046;
signal = dac_o(dataoffset:end);
%signal = dac_data;

%码片搜索 %%采样频率使我们的码片速率的4倍，所以我们的半个码片采样2点  一个码片2046个码元，一个码元采集4次，
samplespercode = round(SampleFreq/(CodeFreq/CodeLength));%每一个PN码循环周期的采样点数，8.184M/(2.046M/2046)
codewithhalf = round((SampleFreq/CodeFreq)/2); %半个码片采样点的数量
codeSerchRange = (samplespercode/codewithhalf);%按照0.5个码片进行搜索，搜索范围 =2
%产生pncode
%%%%前码：pn_code_re 相当于就是我们总的8184个码元，然后如果0到8184，码是0到8183，由于是0.5个码片搜索。所以
%0到4091，个码元，最后一个
%PN码是2046个，就相当于是0到2045，然后前一个码，就是相当于0码元这个位置跑到了2045这码元的位置，他是循环的
%由于我们的码元
%前一个码元就是2045跑到了0码元的位置，然后拼接上后面的码元
%
[pn_code,pn_code_freq] = pn_code_gen(PRN_num1,samplespercode,SampleFreq,CodeFreq,CodeLength); %PRN6
pn_pre = [pn_code_freq((samplespercode - codewithhalf + 1): samplespercode),pn_code_freq(1:(samplespercode - codewithhalf))];
%pn_pre = [pn_code_freq((samplesPerCodecodewidthhalf+1):samplesPerCode),pn_code_freq(1:(samplesPerCodecodewidthhalf))];
pn_cr = pn_code_freq;%当前时刻的pn码相位
%pn_post = [pn_code_freq((codewithhalf + 1):end),pn_code_freq(1,codewithhalf)];
pn_post = [pn_code_freq((codewithhalf+1):end),pn_code_freq(1:codewithhalf)]; %晚码


squ_sum_pre = zeros(1,codeSerchRange); %
squ_sum_cr = zeros(1,codeSerchRange);%
squ_sum_post = zeros(1,codeSerchRange);

sinCarr = zeros(1,samplespercode); %%8184
cosCarr = zeros(1,samplespercode);

temp_offset = zeros(1,samplespercode);
%产生8184点载波
wr_1 = 0;
romaddr_1 = 0;

for dataIndex = 1:samplespercode
 
    wr_1 = wr_1 + CWFreq_word;
    if(wr_1 > 2^32)
        wr_1 = wr_1 - 2^32;
    end
    romaddr_1 = round((wr_1/2^22));
    if(romaddr_1 == 0)
        romaddr_1 = 1; 
    end
   sinCarr(dataIndex) = s_rom(romaddr_1);%%%生成了载波Q路
   cosCarr(dataIndex) = c_rom(romaddr_1);%%生成了载波I路
end

s_t = signal(1:samplespercode); % 取得的第一次8184个数据
signaloffset = 0;

for codeindex  = 1: codeSerchRange
    %1ms的积分
    i_data = s_t.*cosCarr;%i路的数据
    q_data = s_t.*sinCarr;%q路的数据
    temp_pre_i = sum(i_data .* pn_pre);
    temp_cr_i = sum(i_data .* pn_cr);
    temp_post_i = sum(i_data .* pn_post);

    temp_pre_q = sum(q_data .* pn_pre);
    temp_cr_q = sum(q_data .* pn_cr);
    temp_post_q = sum(q_data .* pn_post);

    %功率,I/Q;路
    squ_sum_pre(codeindex) = temp_pre_i^2 + temp_pre_q^2; 
    squ_sum_cr(codeindex) = temp_cr_i^2 + temp_cr_q^2; 
    squ_sum_post(codeindex) = temp_post_i^2 + temp_post_q^2; 

    %%%%数据偏移代替码片偏移，更新源数据，我们1:8184，为什么还要加codewithhalf？每次都是一个码片啊
    signaloffset = signaloffset + samplespercode + codewithhalf  ; %不懂
   % temp_offset(dataIndex) = signaloffset;
    s_t = signal(signaloffset + 1:signaloffset + samplespercode);

end
figure(1);
plot(squ_sum_pre,'r');
hold on;
plot(squ_sum_cr,'b');
hold on ;
plot(squ_sum_post,'g');

1
2
3
4
5
6
7
8
for freqSearchIndex = 1: FreqSearchRange
        freqBin = IF -FreqSearchBand/2 + freqSearchIndex*500;
        freqWord = freqBin*2^32/SampleFreq;
    %产生8184点的载波
        wr=0;
        romaddr=0;
        for dataIndex=1:samplesPerCode
            wr = wr + freqWord;