//Simulateur de processus de Rayleigh //Le code original est un code Matlab qui a été fourni //par J.P Cances (ENSIL) dans le cours de DEA en Traitement du //signal numérique. //Le code délivre un processus aléatoire de Rayleigh dont la Densité //Spectrale de Puissance est déterminée par la vitesse d'un véhicule //en mouvement (spectre de Jakos) //RAZ des variables //clear; //polyfit est une fonction qui permet de déterminer les coefficients //d'un pôlynôme donnée par un vecteur P(X). A l'origine c'est une fonction //matlab qui a été traduite en scilab (trouvé sur le site des //Newsgroup de scilab). //getf('./polyfit.sci'); /////////////////////////////////// //Paramètres d'entrée de la simulation //////////////////////////////////// fd = 0.02; //définition de la fréquence Doppler normalisée fs = 1; //definition de la fréquence d'échantillonnage Ns=1024; //définition du nombre d'échantillon de Rayleigh à produire //Calcul du nombre d'échantillons aléaotoires en entrée N = 8; while (N) if (N<2*fd*Ns/fs) N = 2*N; else break; end end //Nombre de points de la fft inverse (pour lissage) N_inv = ceil(N*fs/(2*fd)); //calcul de l'intervalle fréquentiel après la fft delta_f = 2*fd/N; //calcul de l'intervalle temporel des échantillons de sorties delta_T_inv = 1/fs; //Calcul de deux processus aléatoire gaussien (WGN) I_input_time = rand(1,N,'normal'); Q_input_time = rand(1,N,'normal'); //Calcul des fft des deux processus I_input_freq = fft(I_input_time,1); Q_input_freq = fft(Q_input_time,1); ///////////////////////////////////////////////////// //Calcul de la fonction de transfert du filtre Doppler ///////////////////////////////////////////////////// //Composante continue 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 //utilisation de polyfit pour le calcul de la composante à f = fd k=3; p=polyfit( f(N/2-k:N/2), SEZ(N/2-k:N/2), k); //approximation de la composante à partir des coef donné par polyfit l=size(p,'c'); polyn=0; x=poly(0,"x"); for i=1:l polyn=polyn+p(i)*x^(i-1); end SEZ(N/2+1)=horner(polyn,f(N/2)+delta_f); //mise en forme et affichage de la fonction de transfert SEZ=matrix(SEZ,1,N); xset("window",1);xset("wdim",300,200);xbasc(1);plot2d(SEZ); //Réalisation du filtrage I_output_freq = I_input_freq .* sqrt(SEZ); Q_output_freq = Q_input_freq .* sqrt(SEZ); //insertion de zéros pour le lissage I_temp = [I_output_freq(1:N/2) zeros(1,N_inv-N) I_output_freq(N/2+1:N)]; Q_temp = [Q_output_freq(1:N/2) zeros(1,N_inv-N) Q_output_freq(N/2+1:N)]; //calcul des fft inverses des processus filtrés I_output_time = fft(I_temp,-1); Q_output_time = fft(Q_temp,-1); //calcul du module de l'amplitude au carré de la réponse temporelle for j=1:N_inv r(j) = sqrt( (abs(I_output_time(j)))^2 + (abs(Q_output_time(j)))^2); end //normalisation de la réponse temporelle rms = sqrt( mean (r.*r) ); r = r(1:Ns)/rms; //affichage du processus de Rayleigh xset("window",2);xset("wdim",300,200);xbasc(2);plot2d(r);