clear all;
close all;
echo on
%------------------系統仿真參數
A=1;    %載波振幅
fc=2;   %載波頻率（Hz）
snr=5;  %信噪比dB
N_sample=8;% 基帶信號中每個碼元的的采樣點數

N=10000;  %  碼元數
Ts=1;   % 碼元寬度
df=0.01%頻率分辨率
B=1/Ts;
f_start=fc-B;
f_cutoff=fc+B;
fs=fc*N_sample%系統采樣頻率，即考慮載波后一個碼元內的采樣點數。
ts=Ts/fs; % 系統采樣間隔
t=0:ts:N*Ts-ts;
Lt=length(t);
%-----------畫出調制信號波形及功率譜
% 產生二進制信源
d=sign(randn(1,N));
dd=sigexpand((d+1)/2,fc*N_sample);
gt=ones(1,fc*N_sample); % NRZ波形
d_NRZ=conv(dd,gt);
d_NRZ1=d_NRZ(1:Lt);
pause%畫出單極性NRZ波形及其功率譜
figure(1)
subplot(221); 
plot(t,d_NRZ1);% 畫出單極性NRZ信號波形
axis([0 10 0 1.2]);
xlabel('t');
ylabel('單極性信號');
subplot(222); 
[d_NRZ1f,d_NRZ1,df1,f]=T2F(d_NRZ1,ts,df,fs);%求出單極性NRZ信號功率譜
plot(f,10*log10(abs(fftshift(d_NRZ1f).^2/length(f))));% 畫出單極性NRZ信號功率譜
 axis([-3*B 3*B -50 0]);
xlabel('f');
ylabel('單極性信號PDF');
pause%畫出雙極性NRZ波形及其功率譜
d_sjx=2*d_NRZ-1;%生成雙極性NRZ信號
d_sjx1=d_sjx(1:Lt);
subplot(223); 
plot(t,d_sjx1);% 畫出雙極性NRZ信號波形
axis([0 10 0 1.2]);
xlabel('t');
ylabel('雙極性信號');
subplot(224); 
[d_sjx1f,d_sjx1,df1,f]=T2F(d_sjx1,ts,df,fs);%求出雙極性NRZ信號功率譜
plot(f,10*log10(abs(fftshift(d_sjx1f).^2/length(f))));% 畫出雙極性NRZ信號功率譜
axis([-3*B 3*B -50 0]);
xlabel('f');
ylabel('雙極性信號PDF');
%-------這個圖為figure(1)--------%
 
%---------畫出數字頻帶信號及其功率譜
% 對數字基帶信號進行2ASK調制 
ht=A*sin(2*pi*fc*t);%載波
s_2ask=d_NRZ(1:Lt).*ht;%生成已調信號2ASK
pause%畫出已調信號2ASK及其功率譜
figure(2)
subplot(221);
plot(t,s_2ask);%畫出2ASK信號
axis([0 10 -1.2 1.2]);
xlabel('t');
ylabel('2ASK');
%求2ASK信號功率譜
[s_2askf,s_2ask,df1,f]=T2F(s_2ask,ts,df,fs);%求出單極性NRZ信號功率譜
subplot(222);
plot(f,10*log10(abs(fftshift(s_2askf).^2/length(f))));% 畫出單極性NRZ信號功率譜
axis([-fc-3*B fc+3*B -50 0]);
xlabel('f');
ylabel('2ASK信號PDF');
% 對數字基帶信號進行2FSK調制 
s_2fsk=A*cos(2*pi*fc*t+2*pi*d_sjx(1:Lt).*t);%生成2FSK信號
 pause%畫出已調信號2FSK及其功率譜
subplot(223)
plot(t,s_2fsk);%畫出2FSK波形
axis([0 10 -1.2 1.2]);
xlabel('t');
ylabel('2FSK');
subplot(224)
[s_2fskf,s_2fsk,df1,f]=T2F(s_2fsk,ts,df,fs);%求2FSK信號功率譜
plot(f,10*log10(abs(fftshift(s_2fskf).^2/length(f))));% 畫出2FSK功率譜
axis([-fc-3*B fc+3*B -50 0]);
xlabel('f');
ylabel('2FSK信號PDF');
%-------這個圖為figure(2)--------%
 
% 對數字基帶信號進行2PSK調制 
s_2psk=d_sjx(1:Lt).*ht;%生成2PSK信號
pause%畫出已調信號2PSK及其功率譜
figure(3)
subplot(321)
plot(t,s_2psk);%畫出2PSK波形
axis([0 10 -1.2 1.2]);
xlabel('t');
ylabel('2PSK'); 
subplot(322)
[s_2pskf,s_2psk,df1,f]=T2F(s_2psk,ts,df,fs);%求2PSK信號功率譜
plot(f,10*log10(abs(fftshift(s_2pskf).^2/length(f))));% 畫出2PSK信號功率譜
axis([-fc-3*B fc+3*B -50 0]);
xlabel('f');
ylabel('2PSK信號PDF');
%-----------將2PSK信號送入信道進行傳輸，先生成信道加性高斯白噪聲噪聲
snr_lin=10^(snr/10); %換算成倍數
signal_power=power_x(s_2psk(1:Lt));%求出信號平均功率
noise_power=(signal_power*fs)/(snr_lin*2*(2*B));%求出噪聲方差（噪聲均值為0）
noise_std=sqrt(noise_power);%求出噪聲均方差
noise=noise_std.*randn(1,Lt);%以噪聲均方差作為幅度產生高斯白噪聲
%----------------將已調信號送入信道
pause%畫出信道中的高斯白噪聲及其功率譜
subplot(323)
plot(t,noise(1:Lt));%畫出噪聲
xlabel('t');
ylabel('信道噪聲');
axis([0 10 -3 3]);
[noisef,noise,df1,f]=T2F(noise,ts,df,fs);%求信道噪聲功率譜
subplot(324)
plot(f,10*log10(abs(fftshift(noisef).^2/length(f))));% 畫出信道噪聲功率譜
axis([-fs/2-2 fs/2+2 -50 0]);
xlabel('f');
ylabel('信道噪聲PDF');
r=s_2psk(1:Lt)+noise(1:Lt);%疊加了噪聲的已調信號，相當于將已調信號送入理想信道
pause%畫出加噪后的已調信號2PSK及其功率譜
subplot(325)
plot(t,r);%畫出加噪后的已調信號2PSK
xlabel('t');
ylabel('加噪2PSK信號');
axis([0 10 -3 3]);
[rf,r,df1,f]=T2F(r,ts,df,fs);%求加噪后的已調信號2PSK功率譜
subplot(326)%畫出加噪后已調信號的功率譜
plot(f,10*log10(abs(fftshift(rf).^2/length(f))));% 畫出已調信號2PSK功率譜
axis([-fc-3*B fc+3*B -50 0]);
xlabel('f');
ylabel('加噪2PSK信號PDF');
%-------這個圖為figure(3)--------%
 
%-----在接收端準備進行解調，先通過帶通濾波器
pause%畫出帶通濾波器
[H,f]=bp_f(length(rf),f_start,f_cutoff,df1,fs,1);%經過理想帶通濾波器
figure(4)
subplot(322)
plot(f,abs(fftshift(H)));% 畫出理想帶通濾波器
axis([-fc-3*B fc+3*B -0.2 1.2]);
xlabel('f');
ylabel('理想BPF');
DEM = H.*rf;      %濾波器輸出的頻譜
[dem]=F2T(DEM,fs);%濾波器的輸出波形 
dem1=dem(1:Lt)
pause%經過理想帶通濾波器后的信號波形及功率譜
subplot(323)%經過理想帶通濾波器后的信號波形
plot(t,dem1)%畫出經過理想帶通濾波器后的信號波形
axis([0 10 -4 4]);
xlabel('t');
ylabel('理想BPF輸出信號');
[demf1,dem1,df1,f]=T2F(dem1,ts,df,fs);%求經過理想帶通濾波器后信號功率譜
subplot(324)
plot(f,10*log10(abs(fftshift(demf1).^2/length(f))));% 畫出經過理想帶通濾波器后信號功率譜
axis([-fc-3*B fc+3*B -50 0]);
xlabel('f');
ylabel('理想BPF輸出信號PDF');
%-----進行相干解調，先和本地載波相乘,即混頻
subplot(325)%畫出同頻同相的本地載波
plot(t,ht);
axis([0 10 -1.2 1.2]);
xlabel('t');
ylabel('本地載波');
subplot(326)% 本地載波頻譜
[htf,ht,df1,f]=T2F(ht,ts,df,fs);
plot(f,fftshift(abs(htf)))% 畫出載波頻譜
axis([-fc-3*B fc+3*B 0 15]);
xlabel('f');
ylabel('本地載波頻譜');
%-------這個圖為figure(4)--------%
 
pause%畫出混頻后的信號及功率譜
figure(5)
der=dem1(1:Lt).*ht(1:Lt);%和本地載波相乘，即混頻
subplot(221)%畫出混頻后的波形
plot(t,der);
axis([0 10 -2 2]);
xlabel('t');
ylabel('混頻后的信號');
[derf,der,df1,f]=T2F(der,ts,df,fs);%求混頻后信號的功率譜
subplot(222)
plot(f,10*log10(abs(fftshift(derf).^2/length(f))));%畫出混頻后的功率譜
axis([-fc-3*B fc+3*B -50 0]);
xlabel('f');
ylabel('混頻后信號的PDF');
%---------------再經過低通濾波器
pause%畫出理想低通濾波器
[LPF,f]=lp_f(length(derf),B,df1,fs,1);%求低通濾波器
subplot(224) % 畫出理想低通濾波器
plot(f,fftshift(abs(LPF))); axis([-fc-3*B fc+3*B -0.2 1.2]);
xlabel('f');
ylabel('理想LPF');
pause%混頻信號經理想低通濾波器后的波形及功率譜
DM = LPF.*derf;  %理想低通濾波器輸出信號頻譜
[dm]=F2T(DM,fs); %理想低通濾波器的輸出波形
%-------這個圖為figure(5)--------%
 
figure(6)
subplot(221)
plot(t,dm(1:Lt));%畫出經過低通濾波器后的解調出的波形
axis([0 10 -1.2 1.2]);
xlabel('t');
ylabel('LPF輸出信號');
subplot(222)
[dmf,dm,df1,f]=T2F(dm,ts,df,fs);%求LPF輸出信號的功率譜
plot(f,10*log10(abs(fftshift(dmf).^2/length(f))));%畫出LPF輸出信號的功率譜
axis([-fc-3*B fc+3*B -50 0]);
xlabel('f');
ylabel('LPF輸出信號功率譜');
%-------最后對LPF輸出信號抽樣判決
panjue=zeros(1,N);%建立存儲判決值的矩陣
%抽樣判決，規則：大于0判1，小于判-1
for i=1:N;
    if dm(fc*N_sample*(i-1)+fc*N_sample/2+1)>0;%抽樣判決時刻
        panjue(i)=1;
    else
        panjue(i)=-1;
    end
end
%----畫出判決出的基帶信號波形，并和調制信號比較
rr=sigexpand(panjue,fc*N_sample);
rrt=ones(1,fc*N_sample); % NRZ波形
huifu_NRZ=conv(rr,rrt);
pause%觀察調制信號和解調信號波形
%-------這個圖為figure(6)--------%
 
subplot(224)
plot(t,d_sjx(1:Lt));%調制信號波形
axis([0 10 -1.2 1.2]);
xlabel('t');
ylabel('調制信號');
subplot(223)
plot(t,huifu_NRZ(1:Lt));%解調信號波形
axis([0 10 -1.2 1.2]);
xlabel('t');
ylabel('解調信號');
%-------------------統計誤碼數
numoferr=sum(abs(panjue-d)/2)/N%計算出錯誤碼元數