MATLAB生成C/A码并计算自相关值

坐标上海西南某高校，一门大作业要求我们实现MATLAB生成C/A码并计算自相关值

题目要求：

1. 用Matlab编写程序，生成一个适合PRN 1-32号GPS卫星中任何一个卫星的C/A Code的通用标准函数。长度为1个周期(1ms)，包含1023个chip。
2. 调用这个通用标准函数，选择其中一颗卫星的C/A Code信号，运算它的自相关值，并生成图形。图形例子见下：
   
   事实上，matlab自带了生成C/A码的函数https://www.mathworks.com/matlabcentral/fileexchange/14670-gps-c-a-code-generator

样例代码：
cacode.m

function g=cacode(sv,fs)
% function G=CACODE(SV,FS) 
% Generates 1023 length C/A Codes for GPS PRNs 1-37
%
% 
% g: nx1023 matrix- with each PRN in each row with symbols 1 and 0
% sv: a row or column vector of the SV's to be generated
% 	 valid entries are 1 to 37
% fs: optional number of samples per chip (defaults to 1), fractional samples allowed, must be 1 or greater.
%
% For multiple samples per chip, function is a zero order hold. 
%
%
% For example to generate the C/A codes for PRN 6 and PRN 12 use:
% g=cacode([6 12]),
% and to generate the C/A codes for PRN 6 and PRN 12 at 5 MHz use
% g=cacode([6 12],5/1.023)
% 
%
% For more information refer to the "GPS SPS Signal Specification"
% http://www.navcen.uscg.gov/pubs/gps/sigspec/default.htm
%
% Dan Boschen 12-30-2007
% boschen@loglin.com
% Revision History
% rev 1.0 Dan Boschen 4-15-2007  Initial Release
%
% rev 1.1 Dan Boschen 7-15-2007  Corrected error with taps for PRN30, should be [2,7] was 
% 					incorrect as [1 7]. Thank you Jadah Zak for finding this.
%
%
% rev 1.2 Dan Boschen 12-26-2007 Fixed column index error when ceil ~ L 
%				Thank you Jared Meadows for finding this.
%
% rev 1.3 Dan Boschen 12-30-2007 Changed comment "first order hold" to
% "zero order hold".
%
% rev 1.4 Dan Boschen 6-1-2010   Updated email address in comments
if nargin<2
	fs=1;
end
if (max(sv)>37) || (min(sv)<1) || (min(size(sv))~=1)
	error('sv must be a row or column vector with integers between 1 and 37\n')
end
if fs<1
	error('fs must be 1 or greater\n')
end	
% force integers
testint=round(sv)-sv;
if testint ~= 0 
	warning('non-integer value entered for sv, rounding to closest integer\n');
	sv = round(sv);
end
% table of C/A Code Tap Selection (sets delay for G2 generator)
tap=[2 6;
    3 7;
    4 8;
    5 9;
    1 9;
    2 10;
    1 8;
    2 9;
    3 10;
    2 3;
    3 4;
    5 6;
    6 7;
    7 8;
    8 9;
    9 10;
    1 4;
    2 5;
    3 6;
    4 7;
    5 8;
    6 9;
    1 3;
    4 6;
    5 7;
    6 8;
    7 9;
    8 10;
    1 6;
    2 7;
    3 8;
    4 9
    5 10
    4 10
    1 7
    2 8
    4 10];
% G1 LFSR: x^10+x^3+1
s=[0 0 1 0 0 0 0 0 0 1];
n=length(s);
g1=ones(1,n);	%initialization vector for G1
L=2^n-1;
% G2j LFSR: x^10+x^9+x^8+x^6+x^3+x^2+1
t=[0 1 1 0 0 1 0 1 1 1];
q=ones(1,n);	%initialization vector for G2
% generate C/A Code sequences:
tap_sel=tap(sv,:);
for inc=1:L
    g2(:,inc)=mod(sum(q(tap_sel),2),2);
    g(:,inc)=mod(g1(n)+g2(:,inc),2);
   g1=[mod(sum(g1.*s),2) g1(1:n-1)];
   q=[mod(sum(q.*t),2) q(1:n-1)];
end
%upsample to desired rate
if fs~=1
	%fractional upsampling with zero order hold
	index=0;
	for cnt = 1/fs:1/fs:L
		index=index+1;
		if ceil(cnt) > L   %traps a floating point error in index
			gfs(:,index)=g(:,L);
		else
			gfs(:,index)=g(:,ceil(cnt));
		end
	end 
	g=gfs;
end

关于CA码：
ca码和p码的特点和用途是：C/A码和P码都属于伪随机码，事实上GPS从根本上来说是一个基于码分多址的扩频通信系统，其中的伪随机码一方面起到了将导航电文信号频谱加以扩展的作用，提高其抗干扰的能力，同时也降低了卫星发射信号的功率需求。

而另一方面，C/A码和P码提供了一种测量伪距的手段。得益于它们良好的自相关特性，接收机通过对接收信号与内部复制的伪随机码的相关运算，检测自相关函数的峰值，从而确定接收信号的码相位并计算出从卫星到接收机的空间距离的。

ca码和p码用途是：C/A码是开放给民用的，而P码称为精确定位码，原本是为军用设计的。它们之间的区别在于，P码的码率和循环的周期都要远远高于C/A码。
码率越高意味着一个码的宽度也就越小，例如P码的码率是10.23Mcps而C/A码的码率是1.023Mcps，那么一个码片的长度就是将近300米，可见利用P码测距更为精确。

关于自相关函数在matlab的表示方法：https://blog.csdn.net/scuthanman/article/details/5588138

这里有个文章，可以实现ca码生成并绘制自相关函数图。
但是这个文档不能复制。baidu文库有一篇同意的文章，但是要充会员才能复制。以前baidu还没有这么恶心的。。
这里祝福一下百度公司越办越大（微笑脸）

我用ocr工具识别得到的代码为：

number=input('Please input a number between 1 and 32:');
G1=ones(1,10);
G2=ones(1,10);
m=ones(1,3);
Code2= zeros(1,1);
G2Table=[2,3,4,5,1,2,1,2,3,2,3,5,6,7,8,9,1,2,3,4,5,6,1,4,5,6,7,8,1,2,3,4;
6,7,8,9,9,10,8,9,10,3,4,6,7,8,9,10,4,5,6,7,8,9,3,6,7,8,9,10,6,7,8,9];

for i=1:1023
R(1)=mod(G2(G2Table(1,number))+G2(G2Table(2, number)),2);
rd(i)=mod(R(1)+G1(10),2);
newBit1=[mod(G1(3)+G1(10),2)];
G1=[newBit1 G1(1:9)];
newBit2=[mod(G2(2)+G2(3)+G2(6)+G2(8)+G2(9)+G2(10),2)];
G2=[newBit2 G2(1:9)];
end
cacode=repmat(rd,1,2);

for i=1:2046
if cacode(i)==0
    cacode(i)=1;
else
    cacode(i)=-1;
end
end
k=length(cacode);
xk=fft(cacode,k);
rm=real(ifft(conj(xk).*xk))/2046;
rm=ifftshift(rm);
axis([0,2000,-0.2,1]);
plot(rm)
set(gca,'xtick',[0:80:k-1]);
xlabel('Code Phase')
ylabel('Normalized Correlation Fuction')
title('Autocorrelation of C/A code')