clear all;
close all;

fs=8e5;%抽樣頻率
fm=20e3;%基帶頻率
n=2*(6*fs/fm);
final=(1/fs)*(n-1);
fc=2e5; % 載波頻率
t=0:1/fs:(final);
Fn=fs/2;%耐奎斯特頻率

%=========================================
%用正弦波產生方波
%==========================================
twopi_fc_t=2*pi*fm*t; 
A=1;
phi=0;
x = A * cos(twopi_fc_t + phi);

% 方波
am=1;
x(x>0)=am;
x(x<0)=-1;

figure(1)
subplot(321);
plot(t,x);
axis([0 2e-4 -2 2]);
title('基帶信號');
grid on

car=sin(2*pi*fc*t);%載波
ask=x.*car;%載波調制

subplot(322);
plot(t,ask);
axis([0 200e-6 -2 2]);
title('PSK信號');
grid on;

%=====================================================
vn=0.1;
noise=vn*(randn(size(t)));%產生噪音
subplot(323);
plot(t,noise);
grid on;
title('噪音信號');
axis([0 .2e-3 -1 1]);


askn=(ask+noise);%調制后加噪
subplot(324);
plot(t,askn);
axis([0 200e-6 -2 2]);
title('加噪后信號');
grid on;

%======================================================================
%帶通濾波
%======================================================================
fBW=40e3;
f=[0:3e3:4e5];
w=2*pi*f/fs;
z=exp(w*j);
BW=2*pi*fBW/fs;
a=.8547;%BW=2(1-a)/sqrt(a)
p=(j^2*a^2);
gain=.135;
Hz=gain*(z+1).*(z-1)./(z.^2-(p));

subplot(325);
plot(f,abs(Hz));
title('帶通濾波器');
grid on;


Hz(Hz==0)=10^(8);%avoid log(0)
subplot(326);
plot(f,20*log10(abs(Hz)));
grid on;
title('Receiver -3dB Filter Response');
axis([1e5 3e5 -3 1]);

%濾波器系數
a=[1 0 0.7305];%[1 0 p]
b=[0.135 0 -0.135];%gain*[1 0 -1]
faskn=filter(b,a,askn);
figure(2)
subplot(321);
plot(t,faskn);
axis([0 100e-6 -2 2]);
title('通過帶通濾波后輸出');
grid on;


cm=faskn.*car;%解調
subplot(322);
plot(t,cm);
axis([0 100e-6 -2 2]);
grid on;
title('通過相乘器后輸出');

%===================================================================
%低通濾波器
%==================================================================
p=0.72;
gain1=0.14;%gain=(1-p)/2
Hz1=gain1*(z+1)./(z-(p));

subplot(323);
Hz1(Hz1==0)=10^(-8);%avoid log(0)
plot(f,20*log10(abs(Hz1)));
grid on;
title('LPF -3dB response');
axis([0 5e4 -3 1]);

%濾波器系數
a1=[1 -0.72];%(z-(p))
b1=[0.14 0.14];%gain*[1 1]

so=filter(b1,a1,cm);
so=so*10;%add gain
so=so-mean(so);%removes DC component
subplot(324);
plot(t,so);
axis([0 8e-4 -3.5 3.5]);
title('通過低通濾波器后輸出');
grid on;
%======================================================
%Comparator
%======================================================
High=2.5;
Low=-2.5;
vt=0;%設立比較標準
error=0;
len1=length(so);
for ii=1:len1
  if so(ii) >= vt
  Vs(ii)=High;
  else
  Vs(ii)=Low;
  end
end
Vo=Vs;
subplot(325);
plot (t,Vo), title('解調后輸出信號'),
axis([0 2e-4 -5 5])
grid on;
xlabel('時間 (s)'), ylabel('幅度(V)'),