模拟了北斗卫星发送的数据，用来检测代码，通过搜索到卫星.rar

clc; close all; clear all; SampleFreq = 8.184e6;% 采样率 CodeFreq = 2.046e6;%码片速率 CodeLength = 2046;%码长 PRNnum1 = 6;% PRN码序号 DataLen =210;%数据长度 NHLen = 20;%NH码长度 IF =CodeFreq;%中频 %载波ROM DDS Ac=1024; c=0:1/1023:1; c_rom=round(Ac*cos(2*pi*c)); s_rom=round(Ac*sin(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]; NHCode(NHCode==0)=-1; NHCode = repelem(NHCode,1,DataLen);%按照数据长度扩展DataLen倍 ,% BDData = rand(1,DataLen);%产生10个帧数据随机0~1之间数值，需要转换成-1,1 BDData(BDData<0.5) = -1; BDData(BDData>=0.5) = 1; BDData=repelem(BDData,NHLen); %NH码调制过程 ModNH=NHCode.*BDData;%NH码速率 ModNH=repelem(ModNH,(SampleFreq/CodeFreq)*CodeLength);%升采样到和采样率一样频率 SignalLen=(SampleFreq/CodeFreq)*CodeLength*NHLen*DataLen;%信号的总长度按照采样率来产生 [pn1,pn1_freq]= pn_code_gen(PRNnum1,SignalLen,SampleFreq,CodeFreq,CodeLength); ModPN=ModNH.*pn1_freq;%把NH调制之后的数据进行PN调制 %CW频率控制字 CWFreqWord=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 + CWFreqWord; if(wr > 2^32) wr = wr -2^32; end romaddr = round(wr/2^22); if(romaddr == 0) romaddr =1; end end figure; plot(dac_o(1:160));%调制之后的信号 hold on; plot(ModPN(1:160)*512,'r');%调制完PN码的信号 title('蓝色是dac_o 红色是调制后PN'); %%读文件实现信号捕获 file = fopen('TEST6.BIN','rb'); data = fread(file,'int16'); fclose(file); dac_data = data'; %%数据偏移 dataOffset = 2046; signal = dac_o(dataOffset:end); %%data from file %signal = dac_data; FreqSearchBand = 15e3;%搜索频率范围 hz FreqSearchRange = FreqSearchBand/500; PRNSearchRange = 32; %%码片搜索 samplesPerCode = round(SampleFreq/(CodeFreq/CodeLength));%每个PN码循环周期的采样点数8.184M/(2.046M/2046) codewidthhalf = round((SampleFreq/CodeFreq)/2);%半个码片的采样点数量 codeSearchRange = (samplesPerCode/codewidthhalf);%按照0.5chip 进行搜索，搜索范围是 squ_sum_pre = zeros(FreqSearchRange,codeSearchRange); squ_sum_cr = zeros(FreqSearchRange,codeSearchRange); squ_sum_post = zeros(FreqSearchRange,codeSearchRange); fprintf('('); for PRNIndex = 1:PRNSearchRange %gen pn code [pn_code,pn_code_freq]= pn_code_gen(PRNIndex,samplesPerCode,SampleFreq,CodeFreq,CodeLength);%早码 pn_pre = [pn_code_freq((samplesPerCode-codewidthhalf+1):samplesPerCode),pn_code_freq(1:(samplesPerCode-codewidthhalf))]; pn_cr=pn_code_freq;%当前时刻的pn码相位 pn_post = [pn_code_freq((codewidthhalf+1):end),pn_code_freq(1:codewidthhalf)]; %晚码 sinCarr=zeros(1,samplesPerCode); cosCarr=zeros(1,samplesPerCode); 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; if(wr > 2^32) wr = wr -2^32; end romaddr = round(wr/2^22); if(romaddr == 0) romaddr =1; end sinCarr(dataIndex) = s_rom(romaddr); cosCarr(dataIndex) = c_rom(romaddr); end s_t=signal(1:samplesPerCode);%取得第一次输入的8184点数据 signalOffset =0; for codeindex= 1:codeSearchRange i_data = s_t.*cosCarr; q_data = s_t.*sinCarr; %%1ms 积分 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); %% IQ路功率 squ_sum_pre(freqSearchIndex,codeindex) = temp_pre_i^2 + temp_pre_q^2; squ_sum_cr(freqSearchIndex,codeindex) = temp_cr_i^2 + temp_cr_q^2; squ_sum_post(freqSearchIndex,codeindex) = temp_post_i^2 + temp_post_q^2; %%数据偏移代替码片偏移 signalOffset = signalOffset + samplesPerCode +codewidthhalf; s_t=signal(signalOffset+1:signalOffset+samplesPerCode); end end %判断是否有此prn的相关峰 [peakValue,freqBinIndex]=max(max(squ_sum_cr,[],2));%搜索到最大值的频点 [peakValue,codePhase]=max(squ_sum_cr(freqBinIndex,:));%通过最大值的频点在搜索到相位 preRange = codePhase - codewidthhalf; postRange = codePhase + codewidthhalf; secondSum=squ_sum_cr(freqBinIndex,:); if preRange<=0 secondSum(codeSearchRange+preRange:codeSearchRange)=0; secondSum(1:postRange)=0; elseif postRange>=codeSearchRange secondSum(1:1+postRange-codeSearchRange)=0; secondSum(preRange:codeSearchRange)=0; else secondSum(preRange:postRange)=0; end [secondPeakValue,secondCodePhase] =max(secondSum); if peakValue/secondPeakValue >2.5 fprintf('%02d',PRNIndex);%找到pn码的相关峰并打印PRN编号 else fprintf('.');%未找到相关峰进行打印... end end fprintf(')\n'); %%绘制波形 % figure; % mesh(1:codeSearchRange,1:FreqSearchRange,squ_sum_pre); % title('sum pre'); % %hold on; % figure; % mesh(1:codeSearchRange,1:FreqSearchRange,squ_sum_cr); % title('sum cr'); % %hold on; % figure; % mesh(1:codeSearchRange,1:FreqSearchRange,squ_sum_post); % title('sum post');