% File c7_Jakes.m
% Software given here is to accompany the textbook: W.H. Tranter, 
% K.S. Shanmugan, T.S. Rappaport, and K.S. Kosbar, Principles of 
% Communication Systems Simulation with Wireless Applications, 
% Prentice Hall PTR, 2004.
%
% Generate and test the impulse response of the filter.
%
fd = 100;						% maximum doppler
impw = jakes_filter(fd);		% call to Jakes filter
fs = 16*fd; ts = 1/fs; 			% sampling frequency and time
time = [1*ts:ts:128*ts];		% time vector for plot
subplot(2,1,1)
stem(time,impw,'.'); grid 
xlabel('Time'); ylabel('Impulse Response')
%
% Square the fft and check the power transfer function.
%
[h f] = linear_fft(impw,128,ts);    % generate H(f) for filter
subplot(2,1,2)
plot(f,abs(h.*h)); grid;
xlabel('Frequency'); ylabel('PSD')
%
% Put Gaussian noise through and check the output psd.
%
x = randn(1,1024);                      % generate Gaussian input
y = filter(impw,1,x);                   % filter Gaussian input
[output_psd ff] = log_psd(y,1024,ts);   % log of PSD
figure;
subplot(2,1,1)
plot(ff,output_psd); grid;
axis([-500 500 -50 0])
xlabel('Frequency'); ylabel('PSD')
%
% Filter complex noise and look at the envelope fading.
%
z = randn(1,1024)+i*randn(1,1024);      % generate complex noise
zz = filter(impw,1,z);                  % filter complex noise
time = (0.0:ts:1024*ts);                % new time axis
%
% Normalize output and plot envelope.
%
zz = zz/max(max(abs(zz)));              % normalize to one
subplot(2,1,2)
plot(time(161:480),10*log10(abs(zz(161:480)))); grid
axis([0.1 0.3 -20 0])
xlabel('Time'); ylabel('Log Amplitude')
% End of function file.
% 
% % File: Jakes_filter.m
% % Software given here is to accompany the textbook: W.H. Tranter, 
% % K.S. Shanmugan, T.S. Rappaport, and K.S. Kosbar, Principles of 
% % Communication Systems Simulation with Wireless Applications, 
% % Prentice Hall PTR, 2004.
% %
% function [impw] = jakes_filter(fd)
% % FIR implementation of the Jakes filter (128 points)
% n = 512; nn = 2*n;								% nn is FFT block size
% fs = 0:fd/64:fd;								% sampling frequency = 16*fd
% H = zeros(1,n);									% initialize H(f)
% for k=1:(n/8+1)									% psd for k=1:65
%    jpsd(k)=1/((1-((fs(k))/fd)^2+eps)^0.5);
%    if(jpsd(k))>1000
%       jpsd(k)=1000;
%    end
%    H(k)=jpsd(k)^0.5;							% first 65 points of H
% end
% for k=1:n										% generate negative frequencies
%    H(n+k) = H(n+1-k);
% end
% [inv,time] = linear_fft(H,nn,fd/64);		    % inverse FFT
% imp = real(inv(450:577));						% middle 128 points
% impw = imp.*hanning(128)';						% apply hanning window
% energy = sum(impw.^2);							% compute energy
% impw = impw/(energy^0.5);						% normalize
% % End of function file.
% 
% % File: linear_fft.m
% % Software given here is to accompany the textbook: W.H. Tranter, 
% % K.S. Shanmugan, T.S. Rappaport, and K.S. Kosbar, Principles of 
% % Communication Systems Simulation with Wireless Applications, 
% % Prentice Hall PTR, 2004.
% %
% function [fftx,freq] = linear_fft(x,n,ts)
% % This function takes n (must be even) time domain samples (real or complex)
% % and finds the PSD by taking (fft/n)^2. The two sided spectrum is
% % produced by shifting the PSD. The array freq provides the appropriate
% % frequency values for plotting purposes.
% % By taking 10*log10(psd/max(psd)) the psd is  normalized. Values beow 60db 
% % are set equal to -60 dB.
% y = zeros(1,n);
% for k=1:n
% 	freq(k) =(k-1-(n/2))/(n*ts);
%    y(k) = x(k)*((-1.0)^(k+1));
% end;
% fftx = fft(y)/n;
% % End of function file.
% 
% % File: log_psd.m
% % Software given here is to accompany the textbook: W.H. Tranter, 
% % K.S. Shanmugan, T.S. Rappaport, and K.S. Kosbar, Principles of 
% % Communication Systems Simulation with Wireless Applications, 
% % Prentice Hall PTR, 2004.
% %
% function [logpsd,freq,ptotal,pmax] = log_psd(x,n,ts)
% % This function takes the n time domain samples (real or complex)
% % and finds the psd by taking (fft/n)^2. The two sided spectrum is
% % produced by shifting the psd; The array freq provides the appropriate
% % frequency values for plotting purposes.
% % By taking 10*log10(psd/max(psd)) the psd is  normalized; values beow 60db 
% % are set equal to -60db
% %
% % n must be an even number, preferably a power of  2
% %
% y = zeros(1,n);			% initialize y vector
% %
% h = waitbar(0,'For Loop in PSD Calculation');
% for k=1:n
% 	freq(k) =(k-1-(n/2))/(n*ts);
%    y(k) = x(k)*((-1.0)^k);
%    waitbar(k/n)
% end;
% %
% v = fft(y)/n;
% psd = abs(v).^2;
% pmax = max(psd);
% ptotal = sum(psd);
% logpsd = 10*log10(psd/pmax);
% %
% % Truncate negative values at -60 dB
% %
% for k =1:n
% 	if(logpsd(k)<-60.0)
%       logpsd(k) =-60.0;
% 	end
% end
% close(h)
% % End of function file.