% Fast fading with velocity  Example 6.5

N = 100;                     % number of scatterers
a=randn(N,1);              % create Gaussian amplitude coefficients
a=a/max(a);                 % no coefficent can be larger than 1
th=rand(N,1)*2*pi;         % create uniform phase angles
ph=rand(N,1)*2*pi;
%ph=[zeros(N,1)];
fd=149;   % Doppler
%fd=250;
tmax = 10/fd;               % Maximum time
omega=2*pi*fd;      
t=[0:999]*tmax/999;   % generate timeline
Fs=length(t)/tmax;

X=[zeros(1,length(t))];
Y=[zeros(1,length(t))];
for n=1:N                   % generate the sums for X and Y
 X=X+a(n)*cos(omega*cos(th(n))*t+ph(n));
 Y=Y+a(n)*sin(omega*cos(th(n))*t+ph(n));
end
r=sqrt(X.^2+Y.^2);          %  Calculate the Rayleigh envelope
rdb=20*log10(r);            %  Calculate the envelope in dB
figure;
plot(t*1000,rdb,'k')        % plot
xlabel('time (ms)')
ylabel('envelope')
axis([0 65 -30 10])
set(gca,'units','norm','position',[ .1 .1 .8 .8])
X=[X zeros(1,4*length(X))];
Xf=fft(X);
Xf=fftshift(Xf);
Sd=abs(Xf)/max(abs(Xf));
Sd=Sd.^2;
NXf=length(Xf);
f=[-NXf/2+1:NXf/2]*Fs/(NXf);
figure;
plot(f,abs(Xf)/max(abs(Xf)),'k');
axis([ -Fs/50 Fs/50 0 1.2])
xlabel('Frequency (Hz)')
ylabel('S_d(f)')
set(gca,'units','norm','position',[ .1 .1 .8 .8])