% matlab BPSK 調(diào)制與解調(diào)
clc;
close all;
clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  假定：
%     2倍載波頻率采樣的bpsk信號(hào)
%     調(diào)制速率為載波頻率的 N/2m
%   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=100;%number of bits
fs=10;%采樣頻率 Hz
T=1;%bit interval
t=0:1/fs:m*T-1/fs;%time
fc=3;% carrier frequency
N=fs*T*m;
% N=512;
n=1:1:N;
N0=0.5*randn(1,N); %噪聲
h0=zeros(1,N);
%     30階低通濾波器 h0
f = [0 0.3 0.3 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');
figure(1)
plot(w,h0(1,1:N/2));%低通濾波器，帶寬為采樣頻率的0.3倍,采樣頻率為
for i=1:N/2
    h0(1,N-i+1)=h0(1,i);
end;
%%%%%%%%%   隨機(jī)序列
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信號(hào)
bpsk_m=zeros(1,N);
j=1;k=1;
% sac=sin(2*pi*fc*t+0.05*pi);
coa=cos(2*pi*fc*t+0.05*pi);
for i=1:N
    if(j==(N/m+1))
        j=1;
        k=k+1;
    end;    % 0.05*pi 為初始相位，可以任意改變
%     bpsk_m(1,i)=a(1,k)*sac(i)+a(1,k)*coa(i);
    bpsk_m(1,i)=a(1,k)*coa(i);
    j=j+1;
end;
figure(2)
subplot(2,1,1);
plot(bpsk_m);
xlabel('n');ylabel('bpsk');
pow=mean(abs(bpsk_m).^2)
bpsk_m=bpsk_m+N0;% 信號(hào)加噪聲，模擬過信道
subplot(2,1,2);
plot(bpsk_m);
xlabel('n');ylabel('bpsk_noise');


% 接收處理  用正交本振與信號(hào)相乘，變頻
% bpsk_m1=bpsk_m.*sin(2*pi*fc*t);
bpsk_m2=bpsk_m.*cos(2*pi*fc*t);
%濾波
% tempx=fft(bpsk_m1);
% tempx=tempx.*h0;     %低通濾波
% tempx=ifft(tempx);
% real_x=real(tempx);

% tempx=h0.*fft(bpsk_m2);
 tempx=fft(bpsk_m2);
 figure(3)
 subplot(211);
 plot(abs(tempx));
tempx=tempx.*h0;%低通濾波
subplot(212);
plot(abs(tempx));
tempx=ifft(tempx);
real_x1=real(tempx);

figure(4)
subplot(2,1,1);
% plot(real_x1+real_x,'b'); 
plot(real_x1,'b');
axis([1  N -2.5 2.5]);
grid on;
hold on;
% In=real_x1+real_x;      % 可只取一路，這里取了兩路之和
 In=real_x1;
for i=1:N               % 濾波后整形
    if(In(1,i)>0)       % 判決，得到解調(diào)結(jié)果    
        In(1,i)=1;
    else
        In(1,i)=-1;
    end;
end;
plot(In,'r');

an=zeros(1,m);              
for i=1:m
    an(1,i)=In(1,(i-1)*N/m+N/(2*m));
end;
err=0;
for i=1:m
    if a(i)~=an(i)
        err=err+1;
    end
end
Pe=err/m
subplot(2,1,2);  %  比較誤碼
plot(an,'r*');hold on;
axis([1  m -2 2]);
plot(a,'b^');