% function r = rayleigh_Filter_Model( fd, fs, Ns )
% r = rayleigh(fd,fs,N)
%
% A Rayleigh fading simulator based on Clarke's Model
% Creates a Rayleigh random process with PSD determined
% by the vehicle's speed.
%
% INPUTS:
%   fd = doppler frequency
%        set fd = v*cos(theta)/lambda
%        v = velocity (meters per second)
%        lambda = carrier wavelength (meters)
%        theta = angle w/ respect to tangent (radians).
%   fs = sample frequency (Samples per second)
%   Ns = number of samples of the Rayleigh fading
%        process to produce
%
% OUTPUTS:
%   r  = row vector containing Ns samples of the Rayleigh
%        fading process
%
% Author:  Matthew C. Valenti
%          Mobile and Portable Radio Research Group
%          Virginia Tech
%
% For Academic Use Only


% Remark: fs is the sample rate for received baseband singals. Really, it is 1/T , where fd*T is fading rate
%         Ns is total number to compute 
%         N  is the number of samples between frequencies 0 ~ 2*fd  
%
%   |---------------|----------------------........----------------------------|
%   0              2*fd                                                       fs
%   |<----- N ----->|                                                         |     
%   |<----------------------------- Ns --------------------------------------->|
%

fd=500; %最大多普勒頻移
fs=2e5;  %表示抽樣頻率
Ns=fs/20;  %表示抽樣得到的瑞利衰落信道點(diǎn)數(shù)
N = 8;
while(N)
   if (N < 2*fd*Ns/fs)  % N compare to Ns*(2fd/fs)                            
                                                 N = 2*N;          % N =2.^ m 

   else
      break;
   end
end

% number of ifft points (for smoothing)
N_inv = ceil(N*fs/(2*fd)); %ifft的實(shí)際點(diǎn)數(shù)（考慮平滑作用）

% determine the frequency spacings of signal after FFT   信號(hào)FFT之后的頻率間隔
delta_f = 2*fd/N;

% determine time spacing of output samples  
delta_T_inv = 1/fs; %抽樣點(diǎn)之間的時(shí)間間隔

%%%%%%%%%%% Begin Random Input Generation %%%%%%%%%%%%
% fprintf( 'Generating Input\n');
% generate a pair of TIME DOMAIN gaussian i.i.d. r.v.'s 產(chǎn)生具有N點(diǎn)的高斯隨機(jī)變量
I_input_time = randn(1,N);
Q_input_time = randn(1,N);

% take FFT 做FFT變換
I_input_freq = fft(I_input_time);
Q_input_freq = fft(Q_input_time);

%%%  Generate Doppler Filter's Frequency Response 產(chǎn)生具有經(jīng)典功率譜的成形濾波器 %%%
% fprintf( 'Generating Doppler Filter Function\n');
% filter's DC component
% 
SEZ(1) = 1.5/(pi*fd);

% 0 < f < fd
for j=2:N/2
   f(j) = (j-1)*delta_f;
   SEZ(j) = 1.5/(pi*fd*sqrt(1-(f(j)/fd)^2));
   SEZ(N-j+2) = SEZ(j);
end

% use a polynomial fit to get the component at f = fd
% p = polyfit( f(N/2-3:N/2), SEZ(N/2-3:N/2), 3);
% SEZ(N/2+1) = polyval( p, f(N/2)+delta_f );
% k = N/2 - 1;
% f = fd點(diǎn)的值采用多項(xiàng)式擬合的方法求得
k = 3;
p = polyfit( f(N/2-k:N/2), SEZ(N/2-k:N/2), k );
SEZ(N/2+1) = polyval( p, f(N/2)+delta_f );

%%%%%%%% Perform Filtering Operation 高斯隨機(jī)變量通過(guò)成形濾波器%%%%%%%%%%%%%%
% fprintf( 'Computing Output\n' );
% pass the input freq. components through the filter
I_output_freq = I_input_freq .* sqrt(SEZ);
Q_output_freq = Q_input_freq .* sqrt(SEZ);

% take IFFT 做IFFT變換
I_temp = [I_output_freq(1:N/2) zeros(1,N_inv-N) I_output_freq(N/2+1:N)];
%I_output_time = ifft(I_temp);
I_output_time = real(ifft(I_temp));

Q_temp = [Q_output_freq(1:N/2) zeros(1,N_inv-N) Q_output_freq(N/2+1:N)];
%Q_output_time = ifft(Q_temp);
Q_output_time = real(ifft(Q_temp));

% make vector of times (in milliseconds)
% for j = 1:N_inv
%   t(j) = (j-1)*delta_T_inv*1000;
% end

% take magnitude squared of each component and add together
for j=1:N_inv
   %r(j) = sqrt( (abs(I_output_time(j)))^2 + (abs(Q_output_time(j)))^2);
   r(j)=sqrt(I_output_time(j)^2 + Q_output_time(j)^2);
 end

% normalize and compute rms level
rms = sqrt( mean( r.*r ) );

r=(I_output_time+i*Q_output_time)/rms; %歸一化之后的信道沖激響應(yīng)
r=r(1:Ns);

plot((1:Ns),abs(r));
%s=real(r);
%t=imag(r);
%for ii=1:Ns
%r(ii)=sqrt(s(ii)^2 + t(ii)^2);
%end
%r = r(1:Ns)/rms;