% 4PSK
clc;
close all;
clear;

m=64;
N=256;
n=1:1:N;
h0=zeros(1,N);
N0=0.3*rand(1,N)-0.15; %噪聲


%     30階低通濾波器 h0
f = [0 0.28 0.28 1]; w0 = [1 1 0 0];
b = fir2(30,f,w0);
[h,w] = freqz(b,1,N/2); 
h0(1,1:N/2)=abs(h');
for i=1:N/2
    h0(1,N-i+1)=h0(1,i);
end;

a=rand(1,m);
for i=1:m
    if(a(1,i)<0.25)
        a(1,i)=1;
    else
        if(a(1,i)<0.5)
            a(1,i)=3;
        else
            if(a(1,i)<0.75)
                a(1,i)=5;
            else
                a(1,i)=7;
            end;
        end;
    end;
end;

X_4psk=zeros(1,N);

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

X_4psk_m1=X_4psk.*sin(2*pi*0.25*n+0.15*pi);
X_4psk_m2=X_4psk.*sin(2*pi*0.25*n+0.65*pi);

%濾波
tempx=fft(X_4psk_m1);
tempx=tempx.*h0;     %低通濾波
tempx=ifft(tempx);
X_4psk_m1=real(tempx);


tempx=h0.*fft(X_4psk_m2);
tempx=tempx.*h0;     %低通濾波
tempx=ifft(tempx);
X_4psk_m2=real(tempx);


subplot(3,1,1);
plot(X_4psk_m1);
hold on;
axis([1  N -1.5 1.5]);
grid on;

In1=X_4psk_m1;      
for i=1:N               % 濾波后整形
    if(In1(1,i)>0)       % 判決，得到解調(diào)結(jié)果    
        In1(1,i)=1;
    else
        In1(1,i)=-1;
    end;
end;
plot(In1,'r');

subplot(3,1,2);
plot(X_4psk_m2);
hold on;
axis([1  N -1.5 1.5]);
grid on;

In2=X_4psk_m2;      
for i=1:N               % 濾波后整形
    if(In2(1,i)>0)       % 判決，得到解調(diào)結(jié)果    
        In2(1,i)=1;
    else
        In2(1,i)=-1;
    end;
end;
plot(In2,'r');

b=zeros(1,m);
for i=1:m
    if(In1(1,(i-1)*N/m+N/(2*m))==1&In2(1,(i-1)*N/m+N/(2*m))==1)
        b(1,i)=1;
    end;
    if(In1(1,(i-1)*N/m+N/(2*m))==-1&In2(1,(i-1)*N/m+N/(2*m))==1)
        b(1,i)=3;
    end;
    if(In1(1,(i-1)*N/m+N/(2*m))==-1&In2(1,(i-1)*N/m+N/(2*m))==-1)
        b(1,i)=5;
    end;
    if(In1(1,(i-1)*N/m+N/(2*m))==1&In2(1,(i-1)*N/m+N/(2*m))==-1)
        b(1,i)=7;
    end;
end;
    

subplot(3,1,3);
plot(a,'r*');hold on;
plot(b,'b^');
axis([1  m 0 8]);
grid on;