clc;
close all;
clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  假定：
%     5倍載波頻率采樣的bpsk信號
%     調制速率為在波頻率的 N/2m
%   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=64;%調制碼元個數
N=16384;
n=1:1:N;

N0=1*rand(1,N)-.5; %噪聲

%%%%%%%%%   隨機序列
a=rand(1,m);
for i=1:m
    if(a(1,i)>0.5)
        a(1,i)=1;   
   else
        a(1,i)=-1;
  end;
end;
%%% 生成BPSK信號
bpsk_m=zeros(1,N);
jt=1;k=1;
for i=1:N
    if(jt==(N/m+1))
        jt=1;
        k=k+1;
    end;    % 0.05*pi 為初始相位，可以任意改變
    bpsk_m(1,i)=1.414*(a(1,k)*sin(2*pi*0.2*i)+a(1,k)*cos(2*pi*0.2*i));
    jt=jt+1;
end;

bpsk_m=bpsk_m+N0;% 信號加噪聲，模擬過信道
X_bpsk=bpsk_m;
%% 生成Qpsk信號
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=rand(1,m);
for i=1:m
    if(a(i)<0.25)
        a(i)=-3;
    else
        if(a(i)<0.5)
            a(i)=-1;
        else
            if(a(i)<0.75)
                a(i)=1;
            else
                a(i)=3;
            end;
        end;
    end;
end;
fai=a(1)*pi/4;

X_4psk=zeros(1,N);

jt=1;k=1;
for i=1:N
    if(jt==(N/m+1))
        jt=1;
        k=k+1;
        fai=a(k)*pi/4;
    end;
    
    X_4psk(i)=1.414*sin(2*pi*0.2*i+fai);
    jt=jt+1;
end;
X_4psk=X_4psk+N0;

%%% 生成8PSK信號
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=rand(1,m);
for i=1:m
    if(a(1,i)<0.125)
        a(1,i)=7;
    else
        if(a(1,i)<0.25)
            a(1,i)=-5;
        else
            if(a(1,i)<0.375)
                a(1,i)=3;
            else
                if(a(1,i)<0.5)
                    a(1,i)=-1;
                else
                    if(a(1,i)<0.625)
                        a(1,i)=1;
                    else
                        if(a(1,i)<0.75)
                            a(1,i)=-3;
                        else
                            if(a(1,i)<0.875)
                                a(1,i)=5;
                            else
                                a(1,i)=-7;
                            end;
                        end;
                    end;
                end;
            end;
        end;
    end;
end;

X_8psk=zeros(1,N);
fai=fai+a(1)*pi/8;
jt=1;k=1;
for i=1:N
    if(jt==(N/m+1))
        jt=1;
        k=k+1;
        fai=a(k)*pi/8;
    end;
    X_8psk(1,i)=1.414*sin(2*pi*0.2*i+fai);
    jt=jt+1;
end;

X_8psk=X_8psk+N0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tempx=X_8psk;A=2755;B=3800;
%A=2250;
%B=4250;

temp_fft=fft(tempx);
figure;
plot(abs(temp_fft));
for i=2:A
    temp_fft(i)=0;
    temp_fft(N-i+2)=0;
end;
for i=B:N/2
    temp_fft(i)=0;
    temp_fft(N-i+2)=0;
end;
    temp_fft(1)=0;
    temp_fft(N/2)=0;
    tempx=real(ifft(temp_fft));
for i=1:N/2
    temp_fft(i)=temp_fft(i)*j;
    temp_fft(i+N/2)=temp_fft(i+N/2)*(-j);
end;
temp_x=real(ifft(temp_fft));

figure;
subplot(2,1,1);
plot(abs(tempx+j*temp_x));
hold on;
plot(tempx,'r');
subplot(2,1,2);
plot(abs(temp_fft));

for i=1:N
    fai0(i)=atan2(tempx(i),temp_x(i));
end;

www=2*pi*0.5*(((A+B)/2)/8192);
fai1=mod(www*n,2*pi);
fai0=fai0+pi;
dt=mod(fai0-fai1,2*pi);

qt=diff(dt);
for i=1:N-2
    if(abs(qt(i))>1.9*pi);
        if(i>1)&(i<N-11)
            qt(i)=(qt(i-1)+qt(i+1))/2;
        else
            qt(i)=qt(i+1);
        end;
    end;
end;            

ddt=exp(qt);
figure;
subplot(2,1,1);
plot(dt,'.');
grid on;
subplot(2,1,2);
plot(ddt,'.');
grid on;

var_ddt=var(diff(qt));
ar1=1.045+var_ddt/18;
br1=(0.965-var_ddt/30);

k=1;
i=100;
while(i<N-100)
    if(ddt(i)>ar1)|(ddt(i)<br1)   
        kt=i;
        st=i;
        while(ddt(kt)>ar1)|(ddt(kt)<br1)
            kt=kt-1;
        end;
        while(ddt(st)>ar1)|(ddt(st)<br1)
            ut=ddt(st);
            st=st+1;          
        end;
        utt=abs(mean(dt(kt-4:kt-2))-mean(dt(st+2:st+4)));
        if utt>0.325
            dearta_fai(k)=utt;      
            k=k+1;
        end;
        i=st+1;
    else
        i=i+1;
    end;
end;


d=0:0.001:2*3.2;
q=zeros(size(d));
for i=1:length(q)
    for j=1:k-1
        if (dearta_fai(j)>(d(i)-0.125))&(dearta_fai(j)<(d(i)+0.125))
            q(i)=q(i)+1;
        end;
    end;
end;



figure;
subplot(2,1,1);
plot(dearta_fai,'.');
subplot(2,1,2);
plot(q);

q(1:400)=zeros(1,400);%

d1=10:10:2*3150;
q1=zeros(size(d1));
for i=1:length(q1)
    q1(i)=sum(q(d1(i)-9:d1(i)+9));
end;
figure;
plot(q1);

charge=sum(q1(235:392))/sum(q1);
if charge>0.8
    disp('為BPSK調制方式！');
else
    charge(1)=sum(q1(117:196))/(1+sum(q1(78:235)));
    charge(2)=sum(q1(274:353))/(1+sum(q1(235:392)));
    charge(3)=sum(q1(441:509))/(1+sum(q1(392:548)));
    k=0;
    for i=1:3
        if charge(i)>0.65;
            k=k+1;
        end;
    end;
    if(k>1)
        disp('為4PSK調制方式！');
    else
        charge(1)=sum(q1(58:90))/(1+sum(q1(1:117)));
        charge(2)=sum(q1(137:177))/(1+sum(q1(117:196)));
        charge(3)=sum(q1(215:255))/(1+sum(q1(196:274)));
        charge(4)=sum(q1(294:334))/(1+sum(q1(274:353)));
        charge(5)=sum(q1(372:412))/(sum(q1(353:441))+1);
        charge(6)=sum(q1(450:490))/(sum(q1(441:509))+1);
        charge(7)=sum(q1(528:568))/(sum(q1(509:600))+1);
        k=0;
        for i=1:7
            if charge(i)>0.53;
                k=k+1;
            end;
        end;
        if(k>3)
            disp('為8PSK調制方式！'); 
        else
            disp('無法識別制方式！');
        end;
    end;
end;