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