close all;
N_path=4;                                         
f=1;                                              
f_c=4;                                            
an=1;                                           
N_point=2000;                                     
B=20;                                            
                                                 
T=1/B                                            
time_1 = (linspace (0, N_path, N_point));        
x_min=-0.05;
x_max=0.05;
y_min=0;
                                     
signal_in = sin(2*pi*f*time_1);               
plot (time_1, signal_in, 'b')                  
syms x;                                        
grid on;                                       
xlabel('時間');ylabel('振幅');
title('多徑衰落信道信號')
hold on
signal_1=zeros(1,N_point);                                               
signal_2=zeros(1,N_point);
signal_3=zeros(1,N_point);
n=1
if  n<N_path                                                         
    a=an*exp(-j*2*pi*f_c*n*T);                                             
    
    c=a*exp(-j*2*pi*f_c*time_1(n));                                    

    gg=int(c*sin(pi*B*(x-n*T))/(pi*B*(x-n*T)),x,0,N_path*T);                    
    g=double(gg);                                                          
    s_1=signal_in/cos(2*pi*f_c*time_1(n));
    signal_1=signal_1+g*[zeros(1,n*N_point/N_path*T),s_1(1:end-n*N_point/N_path*T)];   
    signal_out=signal_1;
    plot(time_1,signal_out,'g')                                            
    pause (0.5)                                                            
    
 n=1+n;
    a=an*exp(-j*2*pi*f_c*n*T);                                             
    
    c=a*exp(-j*2*pi*f_c*time_1(n));                                      

    gg=int(c*sin(pi*B*(x-n*T))/(pi*B*(x-n*T)),x,0,N_path*T);                     
    g=double(gg);                                                          
    s_1=signal_in/cos(2*pi*f_c*time_1(n));
    signal_2=signal_2+g*[zeros(1,n*N_point/N_path*T),s_1(1:end-n*N_point/N_path*T)];   
    signal_out=signal_2;
    plot(time_1,signal_out,'g')                                            
    pause (0.5)                                                            
    
n=n+1;                                                        
    a=an*exp(-j*2*pi*f_c*n*T);                                             
    
    c=a*exp(-j*2*pi*f_c*time_1(n));                                      
                                      
    gg=int(c*sin(pi*B*(x-n*T))/(pi*B*(x-n*T)),x,0,N_path*T);                     
    g=double(gg);                                                          
    s_1=signal_in/cos(2*pi*f_c*time_1(n));
    signal_3=signal_3+g*[zeros(1,n*N_point/N_path*T),s_1(1:end-n*N_point/N_path*T)];   
    signal_out=signal_3;
    plot(time_1,signal_out,'g')                                            
    pause (0.5)                                                            
end

signal_out=signal_1+signal_in+signal_2+signal_3;
plot(time_1,signal_out,'r')                                            
legend('輸入信號','多徑信號1','多徑信號2','多徑信號3','輸出信號');

figure;
xlabel('頻率');ylabel('振幅');
title('輸出信號');
grid on;                                                                   
hold on
f_out=-0.5:1/N_point:0.5-1/N_point;                                        
signal_ff=fft(signal_out);                                                 
                                     
signal_f=fftshift(abs(signal_ff));                                       
y_max=max(signal_f);
plot (f_out,signal_f,'r')                                                 


axis([x_min x_max y_min y_max])                                            
hold off


figure;
xlabel('頻率');ylabel('振幅');
title('輸入信號');
grid on;                                                                   
hold on
f_out=-0.5:1/N_point:0.5-1/N_point;                                        

signal1=fft(signal_in);                                                      

signal2=fftshift(abs(signal1));                                            

y_max=max(signal2);
plot (f_out,signal2,'r')                                                   
axis([x_min x_max y_min y_max])                                            

hold off


figure;
xlabel('頻率');ylabel('振幅');
title('輸入輸出信號頻率差值圖');
grid on;                                                                   
hold on
f_out=-0.5:1/N_point:0.5-1/N_point;                                        

 signal_cd=abs(signal2-signal_f)                                                

                                            

y_max=max(signal_cd);
plot (f_out,signal_cd,'r')                                                   
axis([x_min x_max y_min y_max])                                            

hold off


figure;
xlabel('頻率');ylabel('相位');
title('輸出信號相位圖');
grid on;                                                                   
hold on
f_out=-0.5:1/N_point:0.5-1/N_point;                                        
signal_ff=fft(signal_out);                                                 
                                     
signal_f=fftshift(angle(signal_ff));                                       
y_max=max(signal_f);
plot (f_out,signal_f,'r')                                                 


axis([-0.5 0.5 -pi pi])                                            
hold off

figure;
xlabel('頻率');ylabel('相位');
title('輸入信號相位圖');
grid on;                                                                   
hold on
f_out=-0.5:1/N_point:0.5-1/N_point;                                        
signal_ff=fft(signal_in);                                                 
                                     
signal_f=fftshift(angle(signal_ff));                                       
y_max=max(signal_f);
plot (f_out,signal_f,'r')                                                 


axis([-0.5 0.5 -pi pi])                                            
hold off

figure;
xlabel('頻率');ylabel('相位');
title('輸入輸出信號相位差值圖');
grid on;                                                                   
hold on
f_out=-0.5:1/N_point:0.5-1/N_point;                                        

 signal_cd=fftshift(angle(fft(signal_out-signal_in)))                                                
 
 signal_ff=fft(signal_in);                                                 
                                     
signal_f=fftshift(angle(signal_ff));  


                                            

%y_max=max(signal_cd);
plot (f_out,signal_cd-signal_f,'r')                                                   
axis([-0.5 0.5 -pi pi])                                            

hold off