phase_detrending1.m

相位对准技术

t=0.00001-0.00000000001:0.000000000000001:0.00001+0.00000000002;
f1=800000000-100000;
f2=800000000;
f3=800000000+100000;
theta1target=0;
theta2target=0;
theta3target=0;
r1=0.005*randn;
r2=0.005*randn;
r3=0.005*randn;
x1=cos(2*pi*f1*t);
x2=cos(2*pi*f2*t);
x3=cos(2*pi*f3*t);
y1=cos(2*pi*f1*t+0.005*randn(size(t)));
y2=cos(2*pi*f2*t+0.005*randn(size(t)));
y3=cos(2*pi*f3*t+0.005*randn(size(t)));
plot(t,y1,'-.r')
hold on
plot(t,y2,'--b')
plot(t,y3,':g')
hold off
tm=t(1,20000);
theta1tm=2*pi*f1*tm+r1;
theta2tm=2*pi*f2*tm+r2;
theta3tm=2*pi*f3*tm+r3;
theta1tmrem=rem(theta1tm,2*pi);
theta2tmrem=rem(theta2tm,2*pi);
theta3tmrem=rem(theta3tm,2*pi);
kest=round((((theta1tmrem-theta1target)-(theta2tmrem-theta2target))/(2*pi))*(f2/(f2-f1))-(theta2target-theta2tmrem)/(2*pi));
% kest=round(((theta1tmrem-theta2tmrem)/(2*pi))*8000-(0-theta2tmrem)/(2*pi));
test=(2*pi*kest+theta2target-theta2tmrem)/(2*pi*f2)+tm;
kk=[kest-20:kest+20];                   %求k值
tt=(2*pi*kk+theta2target-theta2tmrem)/(2*pi*f2)+tm;            %求相应的t值
% for i=1:41                                                    % 在t值附近求若干个点
%     for j=1:101
%     ttt(i,j)=tt(i)-0.0000000005+0.00000000001*(j-1);
%         end
% end
% for q=1:41
%     for p=1:101
%         tttt(p+101*(q-1))=ttt(q,p);
%     end
% end
tttt=tt(1):0.000000000001:tt(41);
% for i=1:41                                                    % 在t值附近求若干个点
%     for j=1:125
%     ttt(i,j)=tt(i)-0.00000000001*62+0.00000000001*(j-1);
%         end
% end
% for q=1:41
%     for p=1:125
%         tttt(p+101*(q-1))=ttt(q,p);
%     end
% end

theta1tt=rem(2*pi*f1*tttt+r1,2*pi);                             
theta2tt=rem(2*pi*f2*tttt+r2,2*pi); 
theta3tt=rem(2*pi*f3*tttt+r3,2*pi);

for i=1:length(tttt)
   x(i)=theta1tt(i)<pi;
if x(i)==1
        theta1tt(i)=theta1tt(i);
    else
        theta1tt(i)=2*pi-theta1tt(i);
    end
end
for i=1:length(tttt)
   x(i)=theta2tt(i)<pi;
if x(i)==1
        theta2tt(i)=theta2tt(i);
    else
        theta2tt(i)=2*pi-theta2tt(i);
    end
end
for i=1:length(tttt)
   x(i)=theta3tt(i)<pi;
if x(i)==1
        theta3tt(i)=theta3tt(i);
    else
        theta3tt(i)=2*pi-theta3tt(i);
    end
end
% for i=1:length(tttt)
%     if theta1tt(i)>pi
%         theta1tt=theta1tt;
%     else
%         theta1tt=2*pi-theta1tt;
%     end
% end
% for i=1:length(tttt)
%     if theta2tt(i)>pi
%         theta2tt=theta2tt;
%     else
%         theta2tt=2*pi-theta2tt;
%     end
% end
% for i=1:length(tttt)
%     if theta3tt(i)>pi
%         theta3tt=theta3tt;
%     else
%         theta3tt=2*pi-theta3tt;
%     end
% end

% com=abs(theta1tt-theta1target)>pi;
% for i=1:length(theta1tt);
%     if com(i)==1;
%     theta1tt(i)=theta1tt(i);
% else
%     theta1tt(i)=2*pi-theta1tt(i);
%     end
% end
% com=abs(theta2tt-theta2target)>pi;
% for i=1:length(theta2tt);
%     if com(i)==1;
%     theta2tt(i)=theta2tt(i);
% else
%     theta2tt(i)=2*pi-theta2tt(i);
%     end
% end
% com=abs(theta3tt-theta3target)>pi;
% for i=1:length(theta3tt);
%     if com(i)==1;
%     theta3tt(i)=theta3tt(i);
% else
%     theta3tt(i)=2*pi-theta3tt(i);
%     end
% end
E=(theta1tt-theta1target).*(theta1tt-theta1target)+(theta2tt-theta2target).*(theta2tt-theta2target)+(theta3tt-theta3target).*(theta3tt-theta3target);
figure;
plot(tttt,E);