www.gusucode.com > GPS仿真Matlab编程源码程序 > GPS仿真Matlab编程源码程序/Wrelax2.m
function [cPhase,cFreq]= Wrelax2(p, settings,loopCnt)
% function acqResults= Wrelax(p, settings)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clear all
% clc;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% addpath include % The software receiver functions
% addpath geoFunctions % Position calculation related functions
addpath include1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 算法初值设定 =====================================================
M = 1; %% 四阵元
N=settings.samplingFreq*10e-3;%0.05
N1=settings.samplingFreq*11e-3;
Nunit=settings.samplingFreq*1e-3;
samplesPerCode = round(settings.samplingFreq / (settings.codeFreqBasis / settings.codeLength));
% flag=-1;
% count=2*M*N; % 处理4通道20ms数据
timebegin=0;
timeend=N/settings.samplingFreq;%0.05
% m1=zeros(32,1);
% w=zeros(M,32);
% xreal=zeros(24e4,4);
% ximag=zeros(24e4,4);
% data=zeros(1,10);
% data1=zeros(1,10);
% delay1=zeros(10,32);
% cfreq=zeros(10,32);
% peakmatirx=zeros(10,32);
% delta=10;
% L=2^21;
% L=2^(nextpow2(N)+8);
% L1=2^nextpow2(N1);
% cnt=N/N1;
% uniqFftPt = ceil((L + 1) / cnt);
% FftPts = ceil((L1 + 1) / cnt);
t=0:1:N-1;
M1=N/(settings.samplingFreq*1e-3);
L=8*M1*2^nextpow2(N);%15;
% L=2^nextpow2(N);
% m=13;
% delayor=0;
count=0;
snr=-20; % dB
As=10^(snr/20);%1;
% s=zeros(10,10);
% y=zeros(N,2);
%% 仿真数据读取 =====================================================
% fid=fopen( 'realD10.dat','r');
% input_buffer = fread(fid, 2*M*N, 'single');
% fclose(fid);
% fid = fopen('imagD10.dat','r');
% input_buffer2= fread(fid, 2*M*N, 'single');
% fclose(fid);
% for p=1:10
file1=['h:\real_' int2str(p) '.dat'];
fid=fopen(file1,'r');
fseek(fid,4*settings.samplingFreq*loopCnt*1e-3,'bof');
input_buffer = fread(fid, M*N1, 'single');
fclose(fid);
file2=['h:\imag_' int2str(p) '.dat'];
fid = fopen(file2,'r');
fseek(fid,4*settings.samplingFreq*loopCnt*1e-3,'bof');
input_buffer2= fread(fid, M*N1, 'single');
fclose(fid);
xreal=input_buffer(1:N1);
ximag=input_buffer2(1:N1);
% for i=1:M
% xreal(:,i)=input_buffer(1+(i-1)*N1:i*N1);
% ximag(:,i)=input_buffer2(1+(i-1)*N1:i*N1);
% end
% 多天线
% for i=1:M
% xreal(N+1:2*N,i)=input_buffer(1+(i-1+4)*N:(i+4)*N);
% ximag(N+1:2*N,i)=input_buffer2(1+(i-1+4)*N:(i+4)*N);
% end
%% 数据预处理 =====================================================
% snr=-18; % dB
% As=10^(snr/20);%1;
% s=As*signalgen(13,timebegin,settings.samplingFreq,timeend,0.7145e-3,1.6890825e6,codeFreq);%
% x=s.';
x1=xreal+1i*ximag;
x=x1(1:N);
zz=fft(x);
% zz=x(1:6000,1);
% y1=round(real(x.'));
ma=max(abs(real(x1.')));
y0=round(real(x1.')*128/ma);
% y2=round(real(x.')*2^24);
xx=0;
% for i=1:6000:N
% xx=xx+x(i:6000+i-1);
% end
% zz=fft(xx);
% plot(20*log10(abs(fft(zz))))
% z=fft(x(1:N))/max(fft(x(1:N)));
% plot(20*log10(abs(fft(z)/max(abs(fft(z))))))%(-N/2:N/2-1)/N*6e6,
%% 最小功率法求迭代初值 =====================================================
for m=13:13
[codePhase,carrFreq,peak]=acquisition1(y0, settings, m);
if peak>settings.acqThreshold
% for snr=-20:2:0% 卫星数
% for N=N1:N1:N1*8% 卫星数
% if N<10*N1
% M1=N/N1;
% else
% M1=9;
% end
count=count+1;
% uniqFftPt = ceil((L + 1) / cnt);
% timeend=1e-3*N/N1;
% t1=0.257411e-3;
% t2=0.257816e-3;
% fd1=1.690514e6;
% fd2=1.690517e6;
% t1=0.714518e-3;
% t2=0.715248e-3;
% fd1=1.6890825e6;
% fd2=1.6890827e6;
% t1=0.641071e-3;
% t2=0.641771e-3;
% fd1=1.689412e6;
% fd2=1.689415e6;
% t1=0.678741e-3;
% t2=0.679241e-3;
% fd1=1.689412e6;
% fd2=1.689415e6;
% ss1=As*signalgen(m,timebegin,settings.samplingFreq,timeend,t1,fd1,codeFreq);
% ss2=0.85*As*signalgen(m,timebegin,settings.samplingFreq,timeend,t2,fd2,codeFreq);
% ss1=As*signalgen(13,timebegin,settings.samplingFreq,timeend,0.7145e-3,1.6890825e6,codeFreq);%
% [codePhase,carrFreq,peak]=acquisition(y1, settings, m);
% for num=1:100
% if peak>settings.acqThreshold
% k=0;
bn=1;
omiga=0;%0.7205*2*pi*settings.samplingFreq/N;
sm=0;
alpha=0;
tao=0;
flag=0;
tag=0;
% randn('state',100+rand(1)*10);
% Noise=randn(1,N)+1i*randn(M,N);
% zz=fft(ss1).';
% zz=fft(ss1+ss2+Noise).';
% clear Noise
% y=zz;%-alpha1.*sm;
C1=0;%taom=0;
mm=0;
n1=1.686e6;%
% n1=4.306e6;
n2=1.696e6;%
% n2=4.313e6;
step=500;
st=n1:step:n2;
A=zeros(1,length(st));
q=zeros(1,length(st));
for fd=n1:step:n2
% fd=1.6890825e6;
mm=mm+1;
s=As*signalgen(m,timebegin,settings.samplingFreq,timeend,0,fd,settings.codeFreqBasis);%
s=fft(s).';
y1=conj(s).*zz;
y2=abs(fft(y1,L));
% ffts=abs(fft(y1,L));
% [A(mm),q(mm)]=max(ffts);
%-----------------综合Rife方法频率细化---------------%
% loop=0;
% deltaf=settings.samplingFreq/L;
% Z=zeros(L-N,1);
% yt=[y1;Z];
% k0=q(mm);k1=q(mm)-1;k2=q(mm)+1;
% K0=ffts(k0);K1=ffts(k1);K2=ffts(k2);
%
% Kh(1)=2/N*(exp(-1j*(0:L-1)*((k0-0.5)/L)*2*pi)*yt).';
% Kh(2)=2/N*(exp(-1j*(0:L-1)*((k0+0.5)/L)*2*pi)*yt).';
% T=sort(abs(Kh));
% if T(2)/T(1)>1.5
% while loop==0
% if Kh(2)>=Kh(1)
% r=1;
% delta1=r*abs(K2)/(abs(K2)+abs(K0));
% else
% r=-1;
% delta1=r*abs(K1)/(abs(K1)+abs(K0));
% end
% fe=(1-(k0-1+delta1)/L)*settings.samplingFreq;
% dist=mod(fe,deltaf);
% if dist>=1/3*deltaf && dist<=2/3*deltaf
% omiga1=fe*2*pi/settings.samplingFreq;
% loop=1;
% else
% if dist<1/3*deltaf
% r=1;
% deltak=r*(1/2-abs(K2)/(abs(K2+K0)));
% else
% r=-1;
% deltak=r*(1/2-abs(K1)/(abs(K1+K0)));
% end
% k0=k0-deltak;
% k1=k1-deltak;
% k2=k2-deltak;
% K0=2/N*(exp(-1j*(0:L-1)*(1-k0/L)*2*pi)*yt).';
% K1=2/N*(exp(-1j*(0:L-1)*(1-k1/L)*2*pi)*yt).';
% K2=2/N*(exp(-1j*(0:L-1)*(1-k2/L)*2*pi)*yt).';
% Kh(1)=2/N*(exp(-1j*(0:L-1)*((k0-0.5)/L)*2*pi)*yt).';
% Kh(2)=2/N*(exp(-1j*(0:L-1)*((k0+0.5)/L)*2*pi)*yt).';
% end
% end
% else
% delta1=(abs(Kh(2))-abs(Kh(1)))/(abs(Kh(2))+abs(Kh(1))+abs(K0));
% omiga1=(1-(k0+delta1)/L)*2*pi;
% end
%-----------------FFT频率细化---------------%
% delta=q(mm)-1:0.001:q(mm)+1;
% Z=zeros(L-N,1);
% yt=[y1;Z];
% for p=1:length(delta)
% % xf(k)=2/N*sum(x.*exp(-1j*2*pi*(0:1023)*delta(k)/settings.samplingFreq));
% xf(p)=2/N*(exp(-1j*(0:L-1)*(p/L)*2*pi)*yt).';
% end
% [E,I]=max(abs(xf));
% delta1=q(mm)-0.5+(I-1)*0.001;
% % delta1=delta1-1+(I-1)*0.01;
% omiga1=(1-delta1/L)*2*pi;
% clear xf
%--------------------------------------------
% zeros(sf,cnt,1);
% for h=1:cnt
% sf(:,h)=s(N1*(h-1)+1:N1*h);
% yf(:,h)=y(N1*(h-1)+1:N1*h);
% temp(:,h)=fftshift(abs(fft(conj(sf(:,h)).*yf(:,h),L1)));
% [e(h),f(h)]=max(temp(1:FftPts-1,h));
% % [e(h),f(h)]=max(temp(:,h));
% end
% q(mm)=sum(f)/cnt;
% A(mm)=max(sum(temp(1:FftPts-1,:),2));
% [A(mm),q(mm)]=max(sum(temp(1:FftPts-1,:),2));
% omiga1=(q(mm)-1-L1/2)/L1*2*pi;
%-----------------FFT fmin-----------------%
[A(mm),q(mm)]=max(y2);%(1:uniqFftPt-1));
omiga2=(1-(q(mm)-1)/L)*2*pi;
Z=zeros(L-N,1);
yt=[y1;Z];
fun=@(omigam) -abs((exp(-1j*(0:L-1)*omigam)*yt).');
[omiga1,fval,exitflag,output]=fminbnd(fun,omiga2-pi/L,omiga2+pi/L);
clear Z yt
% [A(mm),q(mm)]=max(abs(fft(y1,L)));
% Z=zeros(L-N,1);
% yt=[y1;Z];
% fun=@(p) -abs((exp(-1j*(0:L-1)*(p/L)*2*pi)*yt).');
% kn=fminbnd(fun,q(mm)-1,q(mm)+1);
% omiga1=(1-(kn-1)/L)*2*pi;
%-----------------FFT直接估计---------------%
% [A(mm),q(mm)]=max(fftshift(ffts));
% omiga1=(q(mm)-1-L/2)/L*2*pi;
% if omiga1>0
% omiga1=2*pi-omiga1;
% else
% omiga1=-2*pi-omiga1;
% end
%-----------------FFT czt---------------%
% [A(mm),q(mm)]=max(abs(fft(y1,L)));
% omiga1=(1-(q(mm)-1)/L)*2*pi;
% f1 = omiga1-2*pi/N*0.5; f2 = omiga1+2*pi/N*0.5;
% nn = 2^10;
% % w = exp(-1i*2*pi*(f2-f1)/(nn*fs));
% % b = exp(1i*2*pi*f1/fs);
% w = exp(-1i*(f2-f1)/(nn));
% b = exp(1i*f1);
% B = czt(y1,nn,w,b);
% % a=find(diff(sign(diff(B)))<0)+1;%判断
% % plot(d,B,d(a),B(a),'r*')
% [C,k1]=max(abs(B));
% k1=f1+(k1-1)*(f2-f1)/nn;
% omiga1=k1;
C=max(A);
if C>C1
C1=C;
tao(1)=omiga1*N/(2*pi*settings.samplingFreq);
delay(1)=mod(tao(1)*1e3,1);
omiga(1)=settings.samplingFreq*2*pi*delay(1)*1e-3/N;
fdm=n1+(mm-1)*step;
% fdm=fd;
sm=s.*exp(-1i*omiga(1)*t).';
clear s
mmx=q(mm);
end
end
%-----------------------解扩计算时延----------------------%
% sc1=signalgen1(m,timebegin,settings.samplingFreq,M1*1e-3,delay(1)*1e-3,codeFreq);
% sc2=signalgen1(m,timebegin,settings.samplingFreq,M1*1e-3,0,codeFreq);
% fftNumPt = M1*8*(2^(nextpow2(length(sc1))));
% uniqFftPt = ceil((fftNumPt + 1) / M1);
% fft1=fft(sc1,fftNumPt).*fft(sc2,fftNumPt);
% fft2=fft(sc1,fftNumPt).*fft(sc1,fftNumPt);
% ifft1=ifft(fft1);
% ifft2=ifft(fft2);
%
% [T1,n1]=max(ifft1(1:uniqFftPt-1));
% fun=@(p) -abs(exp(1j*(0:fftNumPt-1)*(p/fftNumPt)*2*pi)*fft1.');
% k1=fminbnd(fun,n1-0.5,n1+0.5);
% [T2,n2]=max(ifft2(1:uniqFftPt-1));
% fun=@(p) -abs(exp(1j*(0:fftNumPt-1)*(p/fftNumPt)*2*pi)*fft2.');
% k2=fminbnd(fun,n2-0.5,n2+0.5);
%
% if k2>k1
% delayt=(k2-k1)/settings.samplingFreq;
% else
% delayt=(6000+k2-k1)/settings.samplingFreq;
% end
% delay(1)=mod(delayt*1e3,1);
% omiga(1)=settings.samplingFreq*2*pi*delay(1)*1e-3/N;
alpha(1)=inv(sm'*sm)*sm'*zz;
% omiga1=omiga(1);
g(1)=(norm(zz-alpha(1)*exp(1i*omiga(1)*t).'))^2; %cost function
g(2)=g(1);
s2=As*signalgen(m,timebegin,settings.samplingFreq,timeend,0,fdm,settings.codeFreqBasis);
s2=fft(s2).';
a=exp(-1i*omiga(1)*t).';
y=zz-alpha(1)*s2.*a;
y1=conj(s2).*y;
y2=abs(fft(y1,L));
[B,k]=max(y2);%(1:uniqFftPt-1));
omiga2=(1-(k-1)/L)*2*pi;
Z=zeros(L-N,1);
yt=[y1;Z];
fun=@(omigam) -abs((exp(-1j*(0:L-1)*omigam)*yt).');
[omiga(2),fval,exitflag,output]=fminbnd(fun,omiga2-pi/L,omiga2+pi/L);
tao(2)=omiga(2)*N/(2*pi*settings.samplingFreq);
delay(2)=mod(tao(2)*1e3,1);
omiga(2)=settings.samplingFreq*2*pi*delay(2)*1e-3/N;
a=exp(-1i*omiga(2)*t).';
alpha(2)=(a'.*s2'*y)/norm(s2.*a)^2;
clear Z yt y1 y2 s2
% alpha1=alpha(2);
% omiga1=omiga(2);
% if delay(1)-delay(2)>1e-4
% temp=delay(2);
% delay(2)=delay(1);
% delay(1)=temp;
% temp=omiga(2);
% omiga(2)=omiga(1);
% omiga(1)=temp;
% end
g(2)=(norm(zz-(alpha*exp(1i*omiga'*t)).'))^2;
d1=inf;d2=inf;
q1=inf;q2=inf;
de=inf;
fm(1)=fdm;
fm(2)=fm(1);
while abs((fm(2)-fm(1)))/fm(1)>3e-5||tag==0
fm(1)=fm(2);
tag=1;
%----------------------- 载波频率细化 -----------------------%
% s11=As*signalgen(m,timebegin,settings.samplingFreq,timeend+1e-3,t1,fd1,codeFreq);%delay(1)*1e-3 delay(2)*1e-3 0.01 0.01
% s22=0.85*As*signalgen(m,timebegin,settings.samplingFreq,timeend+1e-3,t2,fd2,codeFreq);
% Noise1=randn(1,N+settings.samplingFreq*1e-3)+1i*randn(M,N+settings.samplingFreq*1e-3);
% so=s11+s22+Noise1;
so=(x1-mean(x1)).';
clear s11 s22 Noise1
% so=so-mean(so);
% caCode = generateCAcode(m);
%
% codeValueIndex = floor((1 / settings.samplingFreq * (1:M1*samplesPerCode)) / ...
% (1/settings.codeFreqBasis));
%
% longCaCode = caCode((rem(codeValueIndex, 1023) + 1));
longCaCode=signalgen1(m,timebegin,settings.samplingFreq,timeend,0,settings.codeFreqBasis);%delay(1)*1e-3
%--- Remove C/A code modulation from the original signal ----------
% (Using detected C/A code phase)
xCarrier = ...
so(round(delay(1)*samplesPerCode):(round(delay(1)*samplesPerCode) + M1*samplesPerCode-1)) ...
.* longCaCode;
% xCarrier2 = ...
% so(round(delay(2)*samplesPerCode):(round(delay(2)*samplesPerCode) + M1*samplesPerCode-1)) ...
% .* longCaCode;
% xCarrier = xCarrier1+xCarrier2;
% xCarrier = so(1:M1*samplesPerCode).* longCaCode;
%--- Find the next highest power of two and increase by 8x --------
% fftNumPts =8*M1*(2^(nextpow2(length(xCarrier))));%
fftNumPts =2^(nextpow2(length(xCarrier)));
%--- Compute the magnitude of the FFT, find maximum and the *32
%associated carrier frequency
fftxc = abs(fft(xCarrier, fftNumPts));
uniqFftPts = ceil((fftNumPts + 1) / 2);
% [fftMax, fftMaxIndex] = max(fftxc(5 : uniqFftPts-5));
[fftMax, fftMaxIndex] = max(fftxc(1 : uniqFftPts-1));
% [fftMax, fftMaxIndex] = max(fftxc);
% Z=zeros(1,fftNumPts-N);
% yt=[xCarrier,Z];
% fun=@(p) -abs(2/N*yt*exp(-1j*(0:fftNumPts-1)*(p/L)*2*pi).');
% fftMaxIndex1=fminbnd(fun,fftMaxIndex-1,fftMaxIndex+1);
fdm=(fftMaxIndex-1)*settings.samplingFreq/fftNumPts;
% % fftFreqBins = (0 : uniqFftPts-1) * settings.samplingFreq/fftNumPts;
% fftFreqBins = (0 : fftNumPts-1) * settings.samplingFreq/fftNumPts;
%
% % %--- Save properties of the detected satellite signal -------------
% fdm = fftFreqBins(fftMaxIndex);
fm(2)=fdm;
while abs((g(2)-g(1)))/g(1)>3e-5||flag==0
g(1)=g(2);
if flag ==0
flag=1;
% for n=2:-1:1% 迭代次数
% z1=fft(zz.*exp(-1i*2*pi*fdm*t).');
s=As*signalgen(m,timebegin,settings.samplingFreq,timeend,0,fdm,settings.codeFreqBasis);%
s=fft(s).';
y1=conj(s).*zz;
else
y=zz-alpha(2)*s.*a;
y1=conj(s).*y;
end
y2=abs(fft(y1,L));
[B,k]=max(y2);%(1:uniqFftPt-1));
omiga2=(1-(k-1)/L)*2*pi;
Z=zeros(L-N,1);
yt=[y1;Z];
fun=@(omigam) -abs((exp(-1j*(0:L-1)*omigam)*yt).');
[omiga(1),fval,exitflag,output]=fminbnd(fun,omiga2-pi/L,omiga2+pi/L);
tao(1)=omiga(1)*N/(2*pi*settings.samplingFreq);
delay(1)=mod(tao(1)*1e3,1);
omiga(1)=settings.samplingFreq*2*pi*delay(1)*1e-3/N;
a=exp(-1i*omiga(1)*t).';
alpha(1)=(a'.*s'*zz)/norm(s.*a)^2;
clear Z yt y1 y2
% s2=As*signalgen(m,timebegin,settings.samplingFreq,timeend,0,fdm,codeFreq);
% s2=fft(s2).';
% a=exp(-1i*omiga(1)*t).';
y=zz-alpha(1)*s.*a;
y1=conj(s).*y;
y2=abs(fft(y1,L));
% plot(20*log10(abs(fft(y))))(1 : FftPts-1)
% ffts=abs(fft(y1,L));
% [B,k]=max(ffts);
%-----------------综合Rife方法频率细化---------------%
% loop=0;
% deltaf=settings.samplingFreq/L;
% Z=zeros(L-N,1);
% yt=[y1;Z];
% k0=k;k1=k-1;k2=k+1;
% K0=ffts(k0);K1=ffts(k1);K2=ffts(k2);
%
% Kh(1)=2/N*(exp(-1j*(0:L-1)*((k0-0.5)/L)*2*pi)*yt).';
% Kh(2)=2/N*(exp(-1j*(0:L-1)*((k0+0.5)/L)*2*pi)*yt).';
% T=sort(abs(Kh));
% if T(2)/T(1)>1.5
% while loop==0
% if Kh(2)>=Kh(1)
% r=1;
% delta1=r*abs(K2)/(abs(K2)+abs(K0));
% else
% r=-1;
% delta1=r*abs(K1)/(abs(K1)+abs(K0));
% end
% fe=(1-(k0-1+delta1)/L)*settings.samplingFreq;
% dist=mod(fe,deltaf);
% if dist>=1/3*deltaf && dist<=2/3*deltaf
% omiga(n)=fe*2*pi/settings.samplingFreq;
% loop=1;
% else
% if dist<1/3*deltaf
% r=1;
% deltak=r*(1/2-abs(K2)/(abs(K2+K0)));
% else
% r=-1;
% deltak=r*(1/2-abs(K1)/(abs(K1+K0)));
% end
% k0=k0-deltak;
% k1=k1-deltak;
% k2=k2-deltak;
% K0=2/N*(exp(-1j*(0:L-1)*(k0/L)*2*pi)*yt).';
% K1=2/N*(exp(-1j*(0:L-1)*(k1/L)*2*pi)*yt).';
% K2=2/N*(exp(-1j*(0:L-1)*(k2/L)*2*pi)*yt).';
% Kh(1)=2/N*(exp(-1j*(0:L-1)*((k0-0.5)/L)*2*pi)*yt).';
% Kh(2)=2/N*(exp(-1j*(0:L-1)*((k0+0.5)/L)*2*pi)*yt).';
% end
% end
% else
% delta1=(abs(Kh(2))-abs(Kh(1)))/(abs(Kh(2))+abs(Kh(1))+abs(K0));
% omiga(n)=(1-(k0+delta1)/L)*2*pi;
% end
%-----------------FFT频率细化---------------%
% delta=k-1:0.001:k+1;
% Z=zeros(L-N,1);
% yt=[y1;Z];
% for p=1:length(delta)
% % xf(k)=2/N*sum(x.*exp(-1j*2*pi*(0:1023)*delta(k)/settings.samplingFreq));
% xf(p)=2/N*(exp(-1j*(0:L-1)*(p/L)*2*pi)*yt).';
% end
% [E,I]=max(abs(xf));
% delta1=k-1+(I-1)*0.001;
% % delta1=delta1-1+(I-1)*0.01;
% omiga(n)=(1-delta1/L)*2*pi;
% clear xf
%------------------------------------------------
% for h=1:cnt
% sf(:,h)=s2(N1*(h-1)+1:N1*h);
% % y(N1*(h-1)+1:N1*h)=zz(N1*(h-1)+1:N1*h)-alpha1*sf(:,h).*a(N1*(h-1)+1:N1*h);
% yf(:,h)=y(N1*(h-1)+1:N1*h);
% temp(:,h)=fftshift(abs(fft(conj(sf(:,h)).*yf(:,h),L1)));
% [e(h),f(h)]=max(temp(1:FftPts-1,h));
% end
% % k=(sum(f)-min(f)-max(f))/(cnt-2);
% [B,k]=max(sum(temp(1:FftPts-1,:),2));
% omiga(n)=(k-1-L1/2)/L1*2*pi;
%-----------------FFT fmin----------------%
[B,k]=max(y2);%(1:uniqFftPt-1));
omiga2=(1-(k-1)/L)*2*pi;
Z=zeros(L-N,1);
yt=[y1;Z];
fun=@(omigam) -abs((exp(-1j*(0:L-1)*omigam)*yt).');
[omiga(2),fval,exitflag,output]=fminbnd(fun,omiga2-pi/L,omiga2+pi/L);
% [B,k]=max(abs(fft(y1,L)));
% Z=zeros(L-N,1);
% yt=[y1;Z];
% fun=@(p) -abs((exp(-1j*(0:L-1)*(p/L)*2*pi)*yt).');
% kn=fminbnd(fun,k-1,k+1);
% omiga(n)=(1-(kn-1)/L)*2*pi;
%-----------------FFT直接估计---------------%
% [B,k]=max(fftshift(abs(fft(conj(s2).*y,L))));
% omiga(n)=(k-1-L/2)/L*2*pi;
% if omiga(n)>0
% omiga(n)=2*pi-omiga(n);
% else
% omiga(n)=-2*pi-omiga(n);
% end
%-----------------FFT czt---------------%
% [E,k]=max(abs(fft(y1,L)));
% omiga(n)=(1-(k-1)/L)*2*pi;
% f1 = omiga(n)-2*pi/N*0.5; f2 = omiga(n)+2*pi/N*0.5;
% nn = 2^10;
% % w = exp(-1i*2*pi*(f2-f1)/(nn*fs));
% % b = exp(1i*2*pi*f1/fs);
% w = exp(-1i*(f2-f1)/(nn));
% b = exp(1i*f1);
% B = czt(y1,nn,w,b);
% % a=find(diff(sign(diff(B)))<0)+1;%判断
% % plot(d,B,d(a),B(a),'r*')
% [C,k1]=max(abs(B));
% k1=f1+(k1-1)*(f2-f1)/nn;
% omiga(n)=k1;
%-----------------计算时延-----------------%
tao(2)=omiga(2)*N/(2*pi*settings.samplingFreq);
delay(2)=mod(tao(2)*1e3,1);
%-----------------------解扩计算时延----------------------%
% sc1=signalgen1(m,timebegin,settings.samplingFreq,M1*1e-3,delay(n)*1e-3,codeFreq);
% sc2=signalgen1(m,timebegin,settings.samplingFreq,M1*1e-3,0,codeFreq);
% fftNumPt = M1*8*(2^(nextpow2(length(sc1))));
% uniqFftPt = ceil((fftNumPt + 1) / M1);
% fft1=fft(sc1,fftNumPt).*fft(sc2,fftNumPt);
% fft2=fft(sc1,fftNumPt).*fft(sc1,fftNumPt);
% ifft1=ifft(fft1);
% ifft2=ifft(fft2);
%
% [T1,n1]=max(ifft1(1:uniqFftPt-1));
% fun=@(p) -abs(exp(1j*(0:fftNumPt-1)*(p/fftNumPt)*2*pi)*fft1.');
% k1=fminbnd(fun,n1-0.5,n1+0.5);
% [T2,n2]=max(ifft2(1:uniqFftPt-1));
% fun=@(p) -abs(exp(1j*(0:fftNumPt-1)*(p/fftNumPt)*2*pi)*fft2.');
% k2=fminbnd(fun,n2-0.5,n2+0.5);
%
% if k2>k1
% delayt=(k2-k1)/settings.samplingFreq;
% else
% delayt=(6000+k2-k1)/settings.samplingFreq;
% end
% delay(n)=mod(delayt*1e3,1);
omiga(2)=settings.samplingFreq*2*pi*delay(2)*1e-3/N;
a=exp(-1i*omiga(2)*t).';
alpha(2)=(a'.*s'*y)/norm(s.*a)^2;
clear Z yt y1 y2
% alpha1=alpha(n);
% omiga1=omiga(n);
% C1=0;
% for fd=n1:step:n2
% mm=mm+1;
% s11=signalgen(m,timebegin,settings.samplingFreq,timeend,delay(1)*1e-3,fd,codeFreq);
% s22=signalgen(m,timebegin,settings.samplingFreq,timeend,delay(2)*1e-3,fd,codeFreq);
% ss=fft(s11+s22).';
% C=max(real(zz'*ss*inv(ss'*ss)*ss'*zz));
% if C>C1
% C1=C;
% fdm=n1+(mm-1)*step;%*0.1*0.1
% clear ss s11 s22
% end
% end
% errortt=(delay-[t1,t2]*1e3)*1e6;
% if abs(d1)<abs(errortt(1))
% delay(1)=d1*1e-6+t1*1e3;
% omiga(1)=settings.samplingFreq*2*pi*delay(1)*1e-3/N;
% else
% d1=errortt(1);
% end
% if abs(d2)<abs(errortt(2))
% delay(2)=d2*1e-6+t2*1e3;
% omiga(2)=settings.samplingFreq*2*pi*delay(2)*1e-3/N;
% else
% d2=errortt(2);
% end
% end
% de=delay;
% if m==13
% if delay(1)-delay(2)>1e-4
% temp=delay(2);
% delay(2)=delay(1);
% delay(1)=temp;
% temp=omiga(2);
% omiga(2)=omiga(1);
% omiga(1)=temp;
% end
% end
% if abs(delay(1)-t1*1e3)>1e-2%6e-4
% delay(1)=delayor(1);
% else
% if abs(delay(2)-t2*1e3)>1e-2%6e-4
% delay(2)=delayor(2);
% else
% delayor(1)=delay(1);
% delayor(2)=delay(2);
% end
% end
g(2)=(norm(zz-(alpha*exp(1i*omiga'*t)).'))^2;
% errorf=abs(fdm-fd1);
% errort=abs((delay-[t1,t2]*1e3)*1e6);
% if d1<errort(1) && d2<errort(2)
% errort(1)=d1;
% errort(2)=d2;
% break;
% end
% d1=errort(1);d2=errort(2);
% errorf
% errort
end
end
cFreq=fdm;
cPhase=delay(1)*Nunit;
% t_correct=delay(1);
% acqResults.carrFreq(m)=fdm;
% acqResults.codePhase(m)=delay(1)*6000;
% 1111
% result1(num,:)=errort;
% result2(num)=errorf;
% end
% aver1=mean(result1,1);
% aver1=repmat(aver1,100,1);
% aver2=mean(result2,2);
% aver2=repmat(aver2,1,100);
% res1=result1-aver1;
% res2=result2-aver2;
% mse1=sqrt(sum(result1.^2)/num);
% mse2=sqrt(sum(result2.^2)/num);
%
% stdev1=std(result1,0,1);
% stdev2=std(result2,0,2);
%
% mse1re(snr+21,:)=mse1;
% mse2re(snr+21)=mse2;
% mse1re(N/N1,:)=mse1;
% mse2re(N/N1)=mse2;
% figure
% plot(result1(:,1))%,'k:*'
% title('时延误差')
% xlabel('\fontsize{12}次数');
% ylabel('\fontsize{12}误差值(ns)');
% figure
% plot(result1(:,2))
% title('时延误差')
% xlabel('\fontsize{12}次数');
% ylabel('\fontsize{12}误差值(ns)');
% figure
% plot(result2)
% title('载波频率误差')
% xlabel('\fontsize{12}次数');
% ylabel('\fontsize{12}误差值(Hz)');
% figure
% plot(res1(:,1))
% % axis([1 10 -1e-6 1e-6])
% figure
% plot(res1(:,2))
% figure
% plot(res2)
% save file3 acqResults
end
end
% file3=['result_f' int2str(p) '.mat'];
% file4=['result_c' int2str(p) '.mat'];
% fid3=fopen(file3,'a');
% fid4=fopen(file4,'a');
% count3=fwrite(fid3,acqResults.carrFreq,'single');
% count4=fwrite(fid4,acqResults.codePhase,'single');
% fclose(fid3);
% fclose(fid4);
% end
% figure
% plot(-20:2:0,mse1re(1:2:end,1))%,'k:*'
% title('时延估计')
% xlabel('\fontsize{12}SNR(dB)');
% ylabel('\fontsize{12}MSE(ns)');
% figure
% plot(-20:2:0,mse1re(1:2:end,2))%,'k:*'
% title('时延估计')
% xlabel('\fontsize{12}SNR(dB)');
% ylabel('\fontsize{12}MSE(ns)');
% figure
% plot(-20:2:0,mse2re(1:2:end))%,'k:*'
% title('载波频率估计')
% xlabel('\fontsize{12}SNR(dB)');
% ylabel('\fontsize{12}MSE(Hz)');
% figure
% plot(mse1re(:,1))%,'k:*'
% title('时延估计')
% xlabel('\fontsize{12}数据长度(ms)');
% ylabel('\fontsize{12}MSE(ns)');
% figure
% plot(mse1re(:,2))%,'k:*'
% title('时延估计')
% xlabel('\fontsize{12}数据长度(ms)');
% ylabel('\fontsize{12}MSE(ns)');
% figure
% plot(mse2re)%,'k:*'
% title('载波频率估计')
% xlabel('\fontsize{12}数据长度(ms)');
% ylabel('\fontsize{12}MSE(ns)');
% save result.mat mse1re mse2re
% save result1.mat mse1re mse2re
%% 保存结果 =====================================================
% % wopt1=sum(wopt,2);
% % ynew=wopt1'*x;
% % ynew=sum(ynew,1);
% % ynew=zeros(5,2*N);
% ynew=wopt'*x;
% % ynew(1,:)=wopt(:,1)'*x;
% % ynew(2,:)=wopt(:,2)'*x;
% % ynew(3,:)=wopt(:,3)'*x;
% % ynew(4,:)=wopt(:,4)'*x;
% % ynew(5,:)=wopt(:,5)'*x;
% % for i=1:5
% % filename = ['e:\w' int2str(i) '.dat'];
% % fid = fopen(filename,'wb+');
% % filename1 = ['e:\ynew' int2str(i) '.dat'];
% % fwrite(fid,round(real(filename1)*2^5), 'int8');
% % fclose(fid);
% % end
% save ynew.mat ynew wopt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% acquisition the satelliate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Acquisition ===================================================
% % Do acquisition if it is not disabled in settings or if the variable
% % acqResults does not exist.length(Inx)
% for i=4:4
% disp ('Starting acquisiting ...');
% %--- Do the acquisition -------------------------------------------
% disp ('Acquiring satellites...');
% result=round(real(ynew(i,:))*2^5); %
% % result=round(real(ynew(1,:)));
% acqResults = acquisition1(result, settings);
% figure(i);
% plotAcquisition(acqResults);
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%