%%%%%%%%%%%%%%%%%%%%%%%%
%%兩路正弦信號的頻分復用%%
%%%%%%%%%%%%%%%%%%%%%%%%
fs=10000;
f1=50;
f2=100;
t=0:1/fs:1;
x1=sin(2*pi*f1*t);
x2=sin(2*pi*f2*t);

f01=500;
f02=1000;
x01=sin(2*pi*f01*t);
x02=sin(2*pi*f02*t);
y=x1.*x01+x2.*x02;

%時域波形
figure(1);
subplot(4,1,1);plot(t(1:600),x1(1:600));xlabel('sec');
               ylabel('x1');title('原始信號1');
subplot(4,1,2);plot(t(1:600),x2(1:600));xlabel('sec');
               ylabel('x2');title('原始信號2');
               
subplot(4,1,3);plot(t(1:90),x01(1:90));xlabel('sec');
               ylabel('x01');title('載波信號1');
subplot(4,1,4);plot(t(1:90),x02(1:90));xlabel('sec');
               ylabel('x02');title('載波信號2');         

%兩路原始信號頻譜                
figure(2);
N=length(x1);        
fx=fft(x1);          
dfx=fs/N;
n=0:N/2;             
f=n*dfx;            
subplot(2,2,1);plot(f,abs(fx(n+1))*2/N);title('原始信號1頻譜'); 
axis([0 250 0 1]);
xlabel('頻率');
ylabel('幅度');

N=length(x2);        
fx=fft(x2);          
dfx=fs/N;
n=0:N/2;             
f=n*dfx;            
subplot(2,2,2);plot(f,abs(fx(n+1))*2/N);title('原始信號2頻譜'); 
axis([0 250 0 1]);
xlabel('頻率');
ylabel('幅度');

%兩路載波信號頻譜
N=length(x01);        
fx=fft(x01);          
dfx=fs/N;
n=0:N/2;             
f=n*dfx;            
subplot(2,2,3);plot(f,abs(fx(n+1))*2/N);title('載波信號1頻譜'); 
axis([0 1500 0 1]);
xlabel('頻率');
ylabel('幅度');

N=length(x02);        
fx=fft(x02);          
dfx=fs/N;
n=0:N/2;             
f=n*dfx;            
subplot(2,2,4);plot(f,abs(fx(n+1))*2/N);title('載波信號2頻譜'); 
axis([0 1500 0 1]);
xlabel('頻率');
ylabel('幅度');


%傳輸信號頻譜
figure(3);
y=x1.*x01+x2.*x02;
N=length(y);        
fx=fft(y);          
dfx=fs/N;
n=0:N/2;             
f=n*dfx; 
plot(f,abs(fx(n+1))*2/N);title('信道中傳輸信號頻譜'); 
axis([0 1500 0 1]);
xlabel('頻率');
ylabel('幅度');


%設計濾波器通帶截止頻率為400－600Ｈz,阻帶截止頻率為350-650Ｈz, 
Wp1=[400  600]/(fs/2);
Ws1=[350  650]/(fs/2); 
%通帶最大衰減為1db，阻帶最小衰減為５d b
Rp=1;Rs=5;                                  
[N1, Wn1] = BUTTORD(Wp1, Ws1, Rp, Rs); %求出濾波器的階數和固有頻率    
[B1,A1] = butter(N1,Wn1,'bandpass');  %返回濾波器的傳輸函數      
y1 = filter(B1,A1,y);                       

N=length(y1);       %畫出經過濾波器后的第一路高頻載波
fx=fft(y1);      
dfx=fs/N;
n=0:N/2;
f=n*dfx;
figure(4);
subplot(2,1,1);plot(f,abs(fx(n+1))*2/N);title('濾波器濾出的第一路載波');
axis([0 1500 0 0.6]);
xlabel('頻率');
ylabel('幅度');


%以同樣的方法設計第二個帶通濾波器
Wp2=[900 1100]/(fs/2);
Ws2=[850 1150]/(fs/2);   
[N2, Wn2] = BUTTORD(Wp2, Ws2, Rp, Rs);
[B2,A2] = butter(N2,Wn2,'bandpass');
y2 = filter(B2,A2,y);

N=length(y2);           %畫出經過濾波器后的第二路高頻載波
fx=fft(y2);      
dfx=fs/N;
n=0:N/2;
f=n*dfx;
subplot(2,1,2);plot(f,abs(fx(n+1))*2/N);title('濾波器濾出的第二路載波');
axis([0 1500 0 0.5]);
ylabel('幅度');
xlabel('頻率');


%相干解調
Z1=y1.*x01;
Wp3=120/(fs/2);
Ws3=140/(fs/2);
[N3,Wn3]=BUTTORD(Wp3,Ws3,Rp,Rs);
[B3,A3]=butter(N3,Wn3,'low');
y3=filter(B3,A3,Z1);

N=length(y3);           
fx=fft(y3);      
dfx=fs/N;
n=0:N/2;
f=n*dfx;
figure(5);
subplot(2,1,1);plot(f,abs(fx(n+1))*2/N);title('相干解調后的原始信號1頻譜');
axis([0 250 0 0.5]);
ylabel('幅度');
xlabel('頻率');


Z2=y2.*x02;
Wp4=140/(fs/2);
Ws4=160/(fs/2);
[N4,Wn4]=BUTTORD(Wp4,Ws4,Rp,Rs);
[B4,A4]=butter(N4,Wn4,'low');
y4=filter(B4,A4,Z2);

N=length(y4);           
fx=fft(y4);      
dfx=fs/N;
n=0:N/2;
f=n*dfx;
subplot(2,1,2);plot(f,abs(fx(n+1))*2/N);title('相干解調后的原始信號2頻譜');
axis([0 250 0 0.5]);
ylabel('幅度');
xlabel('頻率');


%恢復出的兩路原始信號
figure(6);
subplot(4,1,1);plot(t,y3);xlabel('sec');axis([0 0.1 -0.5 0.5]);
               ylabel('x1');  title('恢復出的原始信號1');
subplot(4,1,2);plot(t,y4);xlabel('sec');axis([0 0.1 -0.47 0.47]);
               ylabel('x2');title('恢復出的原始信號2');
  
               
 %加入噪聲                             
SNRpBit=20;                                            %信噪比
M=2;                                                       %進制數
SNR=SNRpBit/log2(M);
n=awgn(y,SNR-10*log10(0.5)-10*log10(N),'measured',[],'dB');              %在已調信號中加入高斯白噪聲]

%設計濾波器通帶截止頻率為400－600Ｈz,阻帶截止頻率為350-650Ｈz, 
Wp1=[400  600]/(fs/2);
Ws1=[350  650]/(fs/2); 
%通帶最大衰減為1db，阻帶最小衰減為５d b
Rp=1;Rs=5;                                  
[N1, Wn1] = BUTTORD(Wp1, Ws1, Rp, Rs);     
[B1,A1] = butter(N1,Wn1,'bandpass');        
y1 = filter(B1,A1,n);                       

%以同樣的方法設計第二個帶通濾波器
Wp2=[900 1100]/(fs/2);
Ws2=[850 1150]/(fs/2);   
[N2, Wn2] = BUTTORD(Wp2, Ws2, Rp, Rs);
[B2,A2] = butter(N2,Wn2,'bandpass');
y2 = filter(B2,A2,n);


%相干解調
Z1=y1.*x01;
Wp3=120/(fs/2);
Ws3=140/(fs/2);
[N3,Wn3]=BUTTORD(Wp3,Ws3,Rp,Rs);
[B3,A3]=butter(N3,Wn3,'low');
y3=filter(B3,A3,Z1);


Z2=y2.*x02;
Wp4=140/(fs/2);
Ws4=160/(fs/2);
[N4,Wn4]=BUTTORD(Wp4,Ws4,Rp,Rs);
[B4,A4]=butter(N4,Wn4,'low');
y4=filter(B4,A4,Z2);


%兩路加入噪聲后恢復出的原始信號
figure(6);
subplot(4,1,3);plot(t,y3);xlabel('sec');axis([0 0.1 -1.5 1.5]);
               ylabel('x1');  title('加噪后恢復出的原始信號1');
subplot(4,1,4);plot(t,y4);xlabel('sec');axis([0 0.1 -1.5 1.5]);
               ylabel('x2');title('加噪后恢復出的原始信號2');