multipath_doppler.m

多普勒频偏的应用方法的相关程序

%Comming from  Multipath and Doppler effects and Models.Polytechnic 
%University.2005 
%wuguangfu 
%20080918 
 
 
%=============================================================== 
% Note that the receiver signal y(t) remains as a time harmonic signal with the same angular 
% frequency ω as the transmitted signal s(t). Thus, no distortion in wave shape has 
% occurred during the transmission of s(t) through a time invariant multipath environment. 
% However, the magnitude of the signal has been modified. The new magnitude 
% is |H(ω)| which is a function of angular frequency ω . 
 
% ---------------------      figure 2      ---------------------%              
clear all; 
% amplitudes of 7 multipath arrivals 
a=[0.6154 0.7919 0.9218 0.7382 0.1763 0.4057 0.9355]; 
% arrival times of 7 multipath arrivals 
t=[0.9169 0.4103 0.8936 0.0579 0.3529 0.8132 0.0099]; 
i=0; % frequency index 
for omega=0:0.05:100; % angular freuencies 
    multipath_arrival=a.*exp(j*omega*t); 
    i=i+1; 
    abs_H(i)=abs(sum(multipath_arrival)); % the i-th transfer function 
end 
 
omega=0:0.05:100; 
figure (2) 
plot(omega, abs_H) 
ylabel('amplitude of transfer function') 
xlabel('angular freuency') 
title('frequency dependent multipath fading') 
%======================================================== 
 
 
%=============================================================== 
% We will first plot the magnitude of |H(d)| against the distance d using the following 
% matlab code. If the frequency f=1GHz, the wave length is λ =c/ f =0.3 m because the 
% wave speed c=3*10^8m/sec. Let ht=10m, hr=2m. 
% ---------------------      figure 4      ---------------------%        
clear all 
ht=10;hr=2; 
c=3e8;f=1e9;lambda=c/f; 
R=-1; 
d=1:0.5:10000; 
d1=sqrt(d.^2+(ht-hr)^2); 
d2=sqrt(d.^2+(ht+hr)^2); 
a1=exp(j*2*pi.*d1/lambda)./d1; 
a2=R*exp(j*2*pi.*d2/lambda)./d2; 
a=abs(a1+a2); 
ld=log10(d);la=log10(a); 
 
figure (4) 
plot(ld,la); 
xlabel('log10(distance)') 
ylabel('log10(magnitude)') 
title( 'two ray model') 
%============================================================== 
 
 
% ============================================================= 
% Secondly, we plot the magnitude of H(f) against the frequency f for four distances 
% d=50m, 300m, 800m and 2000m using the following matlab code: 
% ---------------------      figure 5      ---------------------%        
clear all 
ht=10;hr=2; 
c=3e8;R=-1; f0=1e8; fi= [1:1:1000];fd=5000000;f= f0+fd*fi; lambda=c./f; 
da=[50,300,800,2000]; 
for i=1:length(da) 
    d=da(i); 
    d1=sqrt(d.^2+(ht-hr)^2); 
    d2=sqrt(d.^2+(ht+hr)^2); 
     %the frequency separation of two adjacent deep fades is in each case 
     %is 1/Td. 
    Td=(d2-d1)/c;   
    a1=exp(j*2*pi*d1./lambda)/d1; 
    a2=R*exp(j*2*pi*d2./lambda)/d2; 
    a(i,:)=abs(a1+a2); 
end 
figure (5) 
 
subplot(2,2,1);plot(f,a(1,:));title('d=50m');ylabel('magnitude') 
subplot(2,2,2);plot(f,a(2,:));title('d=300m');ylabel('magnitude') 
subplot(2,2,3);plot(f,a(3,:)); title('d=800m');xlabel('frequency');ylabel('magnitude') 
subplot(2,2,4);plot(f,a(4,:)); title('d=2000m');xlabel('frequency');ylabel('magnitude') 
%============================================================== 
 
 
% ============================================================= 
% We use the following matlab code to generate the time domain view of transmitted 
% signals and received signals for both cases. From Figure 6, we observe that multipath 
% arrivals cause distortion. The larger the delay spread is, the worse the distortion becomes. 
% ---------------------      figure 6      ---------------------% 
clear all; 
an=[1,0.3,-0.8,0.5,-0.4,0.2];tn=[0,1,2,3,4,5;0,0.1,0.2,0.3,0.4,0.5]; 
signal=[0, zeros(1,0),ones(1,501),zeros(1,1000)]; % transmitted signal 
for k=1:2; %for two cases 
    for i=1:6; 
        ray(i,:)=an(i)*[0, zeros(1,(100*tn(k,i))),ones(1,501),zeros(1,(1000-100*tn(k,i)))]; 
    end 
    y(k,:)=sum(ray(:,1:end));%receiver signal. 
end 
t=((1:1:length(y(1,:)))-1)*10^(-2); 
figure (6) 
subplot(2,2,1);plot(t,signal); 
ylabel('transmitted signal s(t)'); 
title('case 1 & case 2') 
axis([ 0 20 -0.5 1.5]) 
 
subplot(2,2,2);plot(t,y(1,:)); 
ylabel('received signal y(t)'); 
title('case 1: large delay spread') 
 
subplot(2,2,4); 
plot(t,y(2,:)); 
xlabel('Time(us)'); 
ylabel('received signal y(t)'); 
title('case 2: small delay spread') 
%============================================================== 
 
 
% ============================================================= 
% We use the following matlab code to generate the frequency domain view of transmitted 
% signals and received signals for both cases. At first, FFT is employed to implement (3) to 
% find the input spectrum. Secondly, (2) is used to compute the channel transfer functions. 
% Finally, (4) is used to compute the output spectrum. 
% ---------------------      figure 7,8      ---------------------% 
clear all; 
s=[ones(1,10),zeros(1,90)]; % transmitted signal 
s_f=fft(s); 
x=s_f([1:50]); 
y=s_f([51:100]); 
signal_f=[y,x]; %input spectrum 
dt=5/10; % each time interval is 0.01 micro sec 
df=1/(100*dt); 
f_s=df*([0:99]-50);% frequecy vector 
an=[1,0.3,-0.8,0.5,-0.4,0.2]; %amplitudes 
f=f_s; 
w=2*pi*f; 
tn_1=[0,1,2,3,4,5]; % arrival times for case 1 
for i=1:6; 
    h1(i,:)=an(i)*exp(-j*w*tn_1(i)); 
end 
 
h_1=sum(h1(:,1:end));%transfer function 
y_1=h_1.*signal_f;%output spectrum 
tn_2=[0,0.1,0.2,0.3,0.4,0.5]; % arrival times for case 2 
for i=1:6; 
    h2(i,:)=an(i)*exp(-j*w*tn_2(i)); 
end 
 
h_2=sum(h2(:,1:end));%transfer function 
y_2=h_2.*signal_f;%output spectrum 
 
figure(7) 
subplot(2,3,1); 
plot(f_s,abs(signal_f)); 
ylabel('magnitude');title('I/P spectrum') 
subplot(2,3,4); 
plot(f_s,angle(signal_f)); 
ylabel('Phase'); 
xlabel('Frequency(MHz)'); 
 
subplot(2,3,2); 
plot(f,abs(h_1)); 
title('channel 1') 
subplot(2,3,5); 
plot(f,angle(h_1)); 
xlabel('Frequency(MHz)'); 
 
subplot(2,3,3); 
plot(f,abs(h_2)); 
title('channel 2') 
subplot(2,3,6); 
plot(f,angle(h_2)); 
xlabel('Frequency(MHz)'); 
 
 
 
figure(8) 
subplot(2,3,1); 
plot(f_s,abs(signal_f)); 
ylabel('magnitude');title('I/P spectrum') 
subplot(2,3,4); 
plot(f_s,angle(signal_f)); 
ylabel('Phase'); 
xlabel('Frequency(MHz)'); 
 
subplot(2,3,2); 
plot(f,abs(y_1)); 
title('O/P spectrum 1') 
subplot(2,3,5); 
plot(f,angle(y_1)); 
xlabel('Frequency(MHz)'); 
 
subplot(2,3,3); 
plot(f,abs(y_2)); 
title('O/P spectrum 2') 
subplot(2,3,6); 
plot(f,angle(y_2)); 
xlabel('Frequency(MHz)'); 
% ============================================================== 
 
 
%============================================================== 
% From Figures 7, the variation rate (with respect to frequency) of the transfer function is 
% proportional to the delay spread. A larger delay spread causes a faster variation rate in the 
% transfer function. This is illustrated by Figure 9 where the absolute value of the transfer 
% functions for four delay spreads are plotted against the frequency. The code of 
% generating this figure is also shown below. For the case with delay spread=0.2 μ sec, a 
% cycle of variation (from a local peak to the next local peak) is of the order of 5 MHz. 
% Similarly, for the case with delay spread=1 μ sec, 5 μ sec, or 10 μ sec , a cycle of 
% 17 variation (from a local peak to the next local peak) is of the order of 1 MHz, 0.2 MHz, or 
% 0.1 MHz, respectively. 
%  
% ---------------------      figure 9      ---------------------% 
% Figure 9 Absolute values of the transfer functions of four delay spreads. 
clear all; 
N=20; %number of rays 
a=rand(1,N); % amplitudes of N multipath arrivals 
tt=rand(1,N);% times of N multipath arrivals 
f=880:0.005:900; 
delay_spread=0.2; 
t=tt*delay_spread; % arrival times of N multipath arrivals, micro sec 
i=0; % frequency index 
for fi=880:0.005:900; % angular freuencies 
    multipath_arrival=a.*exp(j*2*pi*fi*t); 
    i=i+1; 
    abs_H(i)=abs(sum(multipath_arrival)); % the i-th transfer function 
end 
figure (9) 
subplot(2,2,1) 
plot(f, abs_H) 
ylabel('delay_spread=0.2 micro sec') 
xlabel('freuency, MHz') 
delay_spread=1; 
t=tt*delay_spread; % arrival times of N multipath arrivals, micro sec 
i=0; % frequency index 
for fi=880:0.005:900; % angular freuencies 
    multipath_arrival=a.*exp(j*2*pi*fi*t); 
    i=i+1; 
    abs_H(i)=abs(sum(multipath_arrival)); % the i-th transfer function 
end 
subplot(2,2,2) 
plot(f, abs_H) 
ylabel('delay_spread=1 micro sec') 
xlabel('freuency, MHz') 
delay_spread=5; 
t=tt*delay_spread; % arrival times of N multipath arrivals, micro sec 
i=0; % frequency index 
for fi=880:0.005:900; % angular freuencies 
    multipath_arrival=a.*exp(j*2*pi*fi*t); 
    i=i+1; 
    abs_H(i)=abs(sum(multipath_arrival)); % the i-th transfer function 
end 
subplot(2,2,3) 
plot(f, abs_H) 
ylabel('delay_spread=5 micro sec') 
xlabel('freuency, MHz') 
delay_spread=10; 
t=tt*delay_spread; % arrival times of N multipath arrivals, micro sec 
i=0; % frequency index 
for fi=880:0.005:900; % angular freuencies 
    multipath_arrival=a.*exp(j*2*pi*fi*t); 
    i=i+1; 
    abs_H(i)=abs(sum(multipath_arrival)); % the i-th transfer function 
end 
subplot(2,2,4) 
plot(f, abs_H) 
ylabel('delay_spread=10 micro sec') 
xlabel('freuency, MHz') 
 
%============================================================== 
 
 
%============================================================== 
% multipath signals emit from the source at different angles and 
% arrive at the observer at different angles. Thus, Doppler shifts of different arrivals are 
% usually different from one another 
% |Hωt| is a time varying transfer function which is no longer a time harmonic signal. 
% ---------------------      figure 13      ---------------------% 
clear all; 
N=20 ;% number of multipath arrivals 
a=rand(1,N); %amplitude 
tau=rand(1,N); %arrival time 
f_d=1; 
shift=rand(1,N)*2*f_d-f_d; %Doppler shifts