www.gusucode.com > GPS信号捕获程序，信号源为东方联星Newstar210M采集到的中频数字信号（含GUI界面）- > sourcode/acquisitionfinal.m

    %可改变相关时间search.m performs acquisition on collected data on search7.m
%本地载波的幅度变了
clear;
clc;
% clf;
% ***** initial condition *****
disp('Initial condition');
IF=20491635;                %***IF=20Mhz
worktime=1.0;               %***工作时间是2ms
UnitTime=1.0;               %***单位时间是1s
% fs=16367667;                %***采样频率是16MHz,实验频率
    fs=16363280;
ts=(UnitTime/fs);          %***采样时间间隔是1/fs
Numpms=fix(worktime*fs/1000);            %***每ms中有多少个数据点
    if (rem(Numpms,4)~=0),
        Numpmstemp=Numpms+(4-rem(Numpms,4));    %实际的样点数可能并不为4的整数倍（每个样点用2位二进制数进行量化）
    else
        Numpmstemp=Numpms;
    end
datanum=Numpms;
num=[0:Numpms-1];                           %***数据点的数量
fo=IF-1*fs;                  %***输出频率为fo=fi-1*fs=20-1*16=4MHz,没有多普勒频移

% ***** input data file *****
disp('input data file');
% fid=fopen('C:\MATLAB701\work\GPSRawData.dat', 'r');   %**实验文件
    fid=fopen('C:\MATLAB701\work\GPSData.dat', 'r');
fseek(fid,3000*512, 'bof');               %**从某个位置开始向后读点
buffer1=fread(fid,Numpmstemp,'ubit2');%**之所以有buffer1和buffer是因为数据方向不同
buffer1=buffer1';
buffer=zeros(1,Numpms);
for counter=1:Numpmstemp,    
    switch mod(counter,4)
        case 1,
        temp0=buffer1(counter+3);
        buffer(counter)=temp0;
        case 2,
        temp3=buffer1(counter+1);
        buffer(counter)=temp3;
        case 3,
        temp2=buffer1(counter-1);
        buffer(counter)=temp2;
        case 0,
        temp1=buffer1(counter-3);
        buffer(counter)=temp1;
    end    
    switch buffer(counter)
        case 1
            inputdata(counter)=3;
        case 0
            inputdata(counter)=1;
        case 2
            inputdata(counter)=-1;
        case 3
            inputdata(counter)=-3;
    end
end
clear buffer1;
inputdata=inputdata(1:Numpms);%**按样点数截断读入的数据
finputdata=fft(inputdata);              %**得到各频率谱线
% finputdata=finputdata(1:ceil(length(finputdata)/2));    %取一半的频率值
finputdata(1)=0;                             %**消除直流成分

%******产生本地码，C/A码和RF的乘积后进行环形相关
disp('generate local code and do correlation');
for svID=1:31,
    fprintf('. ');
    CAsample=sampleCAcode(svID,fs,worktime);          %**按照fs采样CA码
    for fideta=1:1:21,                       %**fideta是多普勒频移
        fi(fideta)=fo-10*1000+(fideta-1)*1000;
        lc=CAsample.*(exp(j*2*pi*fi(fideta)*ts*num ));%生成本地载波
        lcf=fft(lc);                                %**做本地载波的fft，即本地波的频谱
%         lcf=lcf(1:ceil(length(lc)/2));    %**取本地波频谱的一半的频率
        out(fideta,:)=abs(ifft(conj(lcf).*(finputdata))); %做环形相关
    end
    [maxAmp(svID) freqshiftNo]=max(max(out'));%每颗星的21个偏移频率中产生最大相关值的那个偏移频率值
    [maxAmp(svID) startpointNo]     =max(max(out ));%每颗星的信号的8184 个频率点,在同21个频率做完相关后,能产生最大的相关值的那个起始位置点
    
  
     outnoise=out;
     outnoise(              : ,startpointNo)=[];
     outnoise(freqshiftNo ,:                 )=[];
     meanvalue=mean(mean(outnoise));
     meanpeak=maxAmp/meanvalue;

end
fprintf('\n');
stem (maxAmp, 'DisplayName', 'maxAmp', 'YDataSource', 'maxAmp'); figure(gcf)