www.gusucode.com > GPS仿真Matlab编程源码程序 > GPS仿真Matlab编程源码程序/Wrelax.m
function [carrFreq,codePhase]= Wrelax(p, settings) % 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*3e-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=M1*8*2^nextpow2(N);%15; % 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'); input_buffer = fread(fid, M*N1, 'single'); fclose(fid); file2=['h:\imag_' int2str(p) '.dat']; fid = fopen(file2,'r'); 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.23e6;% % n1=4.306e6; n2=1.27e6;% % 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 =32*M1*(2^(nextpow2(length(xCarrier))));% %--- Compute the magnitude of the FFT, find maximum and the 8* %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 if (delay(1)-delay(2))>0 && (delay(1)-delay(2))<0.01 temp=delay(1); delay(1)=delay(2); delay(2)=temp; end carrFreq=fdm; codePhase=delay(1)*Nunit; % 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%