Ceci est une ancienne révision du document !
Table des matières
introduction au traitement du signal avec scilab
Introduction au traitement du signal avec Scilab
Représentation d’un signal sinusoïdal dans les domaines temporel et fréquentiel.
// les commentaires sont précédés de deux barres
// on peut utiliser indifféremment majuscules et minuscule, mais « s » n’est pas « S »
// la séparation de la partie entière d’un nombre et des décimales se fait par « . » et non « , »
//
clear
nb_pts=16
// si on ne souhaite pas voir le résultat s’afficher, placer un « ; » à la fin de l’instruction
pas=2e-3
//
t=pas*(0 :1 :nb_pts-1)
// t est un vecteur ligne allant de 0 à nb_pts par pas de 1, multiplié par « pas »
// l’incrément étant par défaut unitaire, il peut être omis ici
//
amp=3 ; f=100 ;
s=amp*sin(2*%pi*f*t)
// les constantes prédéfinies (pi, j, e etc) sont précédées de « % »
//s est un vecteur ligne, de même nombre de points que t
//
// initialisation des paramètres d’affichage
xbasc() ; xset(“font size”, 4);
// affichage
plot2d(t,s) // affichage de s en fonction du temps
// commentaires divers
xtitle (“signal en fonction du temps”,“temps”, “amplitude”);
// initialisation de l’affichage
xbasc(); xset(“font size”,4);
// gestion du premier graphe dans la fenêtre, position (0,0), largeur 1, hauteur 1/2
xsetech([0, 0, 1,1/2]);
// affichage
plot2d(t,s)
// commentaires divers
xtitle (“signal en fonction du temps”,“temps”, “amplitude”);
//
// gestion du second graphe dans la fenêtre, position (0,1/2), largeur 1, hauteur 1/2
xsetech([0, 1/2, 1,1/2]);
// affichage
plot2d(s)
// commentaires divers
xtitle (“signal en fonction du rang des échantillons”,“rang des échantillons”, “amplitude”);
clear
//
// définition des constantes, nombre de points, période d’échantillonnage et fréquence du signal
N=1000 ; Te=0.1e-3 ; F=100 ;
//
// description du vecteur temps et vecteur signal
t=Te*(0:N-1);s=5*sin(2*%pi*F*t);
//
// affichage
xbasc() ; xset(“font size”,5);
plot2d(t,s); xtitle(“”,”temps (s)”,”amplitude (V)”)
// définition d'une variable fréquence pour l'affichage
// le premier point correspond au continu
// le dernier à la fréquence d'échantillonnage au pas près
f=1/(N*Te)*( 0 : N-1) ;
//
// calcul de la transformation
sf=fft(s,-1) ;
//
// initialisation de l'affichage
xbasc(); xset(“font size”,4);
//
// affichage du module brut de forme
xsetech([0,0,1,1/3]) ; plot2d(abs(sf)) ;
xtitle(“module brut de la fft en fonction des indices”,“rang des indices”,“amplitude” );
//
// affichage du module corrigé de 0 à fe/2
xsetech([0,1/3,1,1/3]) ; plot2d(f(1:N/2), 1/N*abs(sf(1:N/2))) ;
xtitle(“module corrigé en fonction de la fréquence de 0 à fe/2”,“fréquence (Hz)”,“amplitude (Vs)”);
//
// zoom sur la raie, de la fréquence nulle à fe/5
xsetech([0,2/3,1,1/3]) ; plot2d(f(1:N/20),1/N*abs(sf(1: N/20))) ;
xtitle(“module corrigé et zoomé”,“fréquence (Hz)”,“amplitude (Vs)”);
Représentation d’un signal périodique quelconque
clear
// définition des constantes : nombre d'échantillons par symbole, nombre de périodes
// nombre total de points, fréquence d'échantillonnage
Nech_symb=10 ; Nsymb=12 ; N=Nech_symb*Nsymb ; Fe=1e3 ;
//
// variable temps
t=1/Fe*[0:N-1];
//
// écriture d'un symbole +/- 5
symb=5*[ ones(1, Nech_symb/2), -1*ones(1, Nech_symb/2)];
//
// répétion du symbole pour constituer le signal
S=matrix(symb' *ones(1, Nsymb), 1 , N);
//
// affichage
xbasc(); xset (“font size”, 4);
plot2d2(t,S, rect=[0, -6, .12, 6]) , xtitle(“signal carré”, “temps (s)”, “amplitude”);
Modulation d’amplitude
**__Modulation sans porteuse__**
clear
// définition des constantes
N=500 ; fe=500e3 ; fp=50e3 ; finf=4e3;
//
// description des vecteurs temps, fréquence (pour l’affichage de la fft) et signal
t=(0:N-1)/fe ; f=(0:N-1)/N*fe ;
sinf=2*sin(2*%pi*finf*t); spor=5*sin(2*%pi*fp*t); smod=spor.*sinf;
//
// initialisation de l'affichage
xbasc() ; xset( “font size”, 4);
//
// affichage du signal informatif
xsetech([0,0,1,1/3]) ; plot2d(t,sinf) ; xtitle(“signal informatif”,“temps (s)”,“amplitude (V)”);
//
// affichage du signal porteur
xsetech([0,1/3,1,1/3]) ; plot2d(t,spor) ;
xtitle(“signal porteur”,“temps (s)”,“amplitude (V)”);
//
// affichage du signal modulé
xsetech([0,2/3,1,1/3]) ; plot2d(t, smod) ;
xtitle(“signal modulé”,“temps (s)”, “amplitude (V)”);
// calcul des modules des fft
Sinf=1/N*abs(fft(sinf,-1)) ; Spor=1/N*abs(fft(spor,-1)) ; Smod=1/N*abs(fft(smod,-1)) ;
//
// initialisation de l'affichage
xbasc() ; xset( “font size”, 4);
//
//
//
// affichage du signal informatif
xsetech([0,0,1,1/3]) ; plot2d3(f(1: .12*N),Sinf(1: .12*N), rect=[0, 0, 6e4, 3]) ;
xtitle(“spectre du signal informatif”,“fréquence (Hz)”,“amplitude (Vs)”);
//
// affichage du signal porteur
xsetech([0,1/3,1,1/3]) ; plot2d3(f(1: .12*N),Spor(1: .12*N), rect=[0, 0, 6e4, 3]) ;
xtitle(“spectre du signal porteur”,“fréquence (Hz)”,“amplitude (Vs)”);
//
// affichage du signal modulé
xsetech([0,2/3,1,1/3]) ; plot2d3(f(1: .12*N),Smod(1: .12*N), rect=[0, 0, 6e4, 3]) ;
xtitle(“spectre du signal modulé”,“fréquence (Hz)”,“amplitude (Vs)”);
clear ;
// définition des constantes
N=500 ; fe=500e3 ; fp=50e3 ; finf1=4e3; finf2=2e3; finf3=5e3;
//
// description des vecteurs temps, fréquence (pour l’affichage de la fft) et signal
t=(0:N-1)/fe ; f=(0:N-1)/N*fe ;
//
sinf=2*sin(2*%pi*finf2*t)+1*sin(2*%pi*finf1*t)+1.5*sin(2*%pi*finf3*t);
spor=5*sin(2*%pi*50e3*t);
smod=spor.*sinf;
//
// initialisation de l'affichage
xbasc() ; xset( “font size”, 4);
//
// affichage du signal informatif
xsetech([0,0,1,1/2]) ; plot2d(t,sinf) ; xtitle(“signal informatif”,“temps (s)”,“amplitude (V)”);
//
// affichage du signal modulé
xsetech([0,1/2,1,1/2]) ; plot2d(t, smod) ;
xtitle(“signal modulé”,“temps (s)”, “amplitude (V)”);
// calcul des modules des fft
Smod=1/N*abs(fft(smod,-1)); Sinf=1/N*abs(fft(sinf,-1));
//
// initialisation de l'affichage
xbasc() ; xset( “font size”, 4);
// affichage du signal informatif
xsetech([0,0,1,1/2]) ; plot2d3(f(1: .12*N),Sinf(1: .12*N), rect=[0, 0, 6e4, 1.5]) ;
xtitle(“spectre du signal informatif”,“fréquence (Hz)”,“amplitude (Vs)”);
//
// affichage du signal modulé
xsetech([0,1/2,1,1/2]) ; plot2d3(f(1: .12*N),Smod(1: .12*N)) ;
xtitle(“spectre du signal modulé”,“fréquence (Hz)”,“amplitude (Vs)”);
clear
// définition des constantes
N=500 ; fe=500e3 ; fp=50e3 ; finf1=4e3; finf2=10e3; finf3=15e3;
//
// description du vecteur temps
t=(0:N-1)/fe ;
//
// les fréquences en vecteur colonne
frequence=[4e3; 10e3; 15e3];
//
// les amplitudes en vecteur ligne
amplitude=[1, 2, 1.5] ;
//
// calcul du signal modulé
sinf_2=amplitude*sin(2*%pi*frequence*t);
//
//affichage
xbasc(); plot2d(t,sinf)
**__Modulation en bande latérale unique__**
clear
// définition des constantes
N=5000 ; fe=500e3 ; fp=50e3 ; finf1=4e3; finf2=3e3; finf3=5e3;
//
// description des vecteurs temps, fréquence (pour l’affichage de la fft) et signal
t=(0:N-1)/fe ; f=(0:N-1)/N*fe ;
//
// définition des différents signaux informatif, porteur et double bande
sinf=2*sin(2*%pi*finf2*t)+1*sin(2*%pi*finf1*t)+1.5*sin(2*%pi*finf3*t);
spor=5*sin(2*%pi*50e3*t);
smod=spor.*sinf;
//
// calcul de la fonction de transfert du filtre et des valeurs de cette fonction sur la gamme de fréquence
G=analpf(8,'cheb1',[.1 0],2*%pi*(fp-finf2));
gain=freq( G(2), G(3), %i*2*%pi*f);
//
// filtrage par multiplication dans le domaine fréquentiel
Smod=fft(smod,-1);
Sblu=Smod.*gain;
//
//affichage dans le domaine spectral
xbasc() ; xset( “font size”, 4);
//
xsetech([0,0,1,1/3]) ; plot2d(f(1:N/5),1/N*abs(Smod(1:N/5))) ;
xtitle(“signal double bande”,“fréquence (Hz)”,“amplitude (Vs)”);
//
xsetech([0,1/3,1,1/3]) ; plot2d(f(1:N/5), abs(gain(1:N/5))) ;
xtitle(“filtre en échelles linéaires”,“fréquence (Hz)”,“gain”);
//
//
xsetech([0,2/3,1,1/3]) ; plot2d(f(1:N/5), 1/N*abs(Sblu(1:N/5))) ;
xtitle(“signal blu”,“fréquence (Hz)”,“amplitude (Vs)”);
// définition des paramètres d'affichage de 10 kHz à 100 kHz
Nmin=N/50 ; Nmax=N/5; f_aff=f(Nmin:Nmax);
//
Smod_dB=20*log10(1/N*abs(Smod(Nmin:Nmax)));
gain_dB=20*log10( abs(gain(Nmin:Nmax)));
Sblu_dB= 20*log10(1/N*abs(Sblu(Nmin:Nmax)));
//
// affichage en dB
xbasc(); xset( “font size”, 4);
//
xsetech([0,0,1,1/3]) ; plot2d(f_aff ,Smod_dB, logflag=[“ln”] , rect=[1e4, -320, 1e5, 20] ) ;
xtitle(“signal double bande”,“fréquence (Hz)”,“amplitude (dB)”);
//
xsetech([0,1/3,1,1/3]) ; plot2d(f_aff, gain_dB , logflag=[“ln”] , rect=[1e4, -90, 1e5, 10] ) ;
xtitle(“filtre”,“fréquence (Hz)”,“gain (dB)”);
//
xsetech([0,2/3,1,1/3]) ; plot2d(f_aff, Sblu_dB, logflag=[“ln”] , rect=[1e4, -380, 1e5, 20] ) ;
xtitle(“signal blu”,“fréquence (Hz)”,“amplitude (dB)”);
//
// retour dans l'espace temporel
Sblu=Sblu(1:N/2+1);
Sblu=[Sblu(1 :N/2), conj(Sblu(N/2+1 : -1 : 2))];
sblu=fft(Sblu,1);
//
// affichage dans le domaine temporel
xbasc() ; xset( “font size”, 4);
xsetech([0,0,1,1/3]) ; plot2d(t (1:N/5) ,sinf(1:N/5) ) ;
xtitle(“signal informatif”,“temps (s)”,“amplitude (V)”);
//
xsetech([0,1/3,1,1/3]) ; plot2d(t(1:N/5) ,smod(1:N/5) ) ;
xtitle(“signal double bande”,“temps (s)”, “amplitude (V)”);
//
xsetech([0,2/3,1,1/3]) ; plot2d(t(1:N/5) ,sblu(1:N/5) ) ;
xtitle(“signal blu”,“temps (s)”, “amplitude (V)”);
// définition des constantes
N=5000 ; fe=500e3 ; fp=50e3 ; finf1=4e3; finf2=3e3; finf3=5e3;
//
// description du vecteur temps
t=(0:N-1)/fe ;
//
// définition des différents signaux informatif, modulé bande supérieure et inférieure
s_inf=2*sin(2*%pi*(finf2)*t)+1*sin(2*%pi*(finf1)*t)+1.5*sin(2*%pi*(finf3)*t);
sblu_inf=5*cos(2*%pi*(fp-finf2)*t)+2.5*cos(2*%pi*(fp-finf1)*t)+3.75*cos(2*%pi*(fp-finf3)*t);
sblu_sup=-5*cos(2*%pi*(fp+finf2)*t)-2.5*cos(2*%pi*(fp+finf1)*t)-3.75*cos(2*%pi*(fp+finf3)*t);
s_dif=sblu_sup-sblu_inf;
//
xbasc() ; xset( “font size”, 4);
xsetech([0,0,1,1/4]) ; plot2d(t (1:N/5) ,s_inf(1:N/5) ) ;
xtitle(“signal informatif”,“temps (s)”,“amplitude (V)”);
//
xsetech([0,1/4,1,1/4]) ; plot2d(t(1:N/5) ,sblu_inf(1:N/5) ) ;
xtitle(“signal blu bande inférieure”,“temps (s)”, “amplitude (V)”);
//
xsetech([0,2/4,1,1/4]) ; plot2d(t(1:N/5) ,sblu_sup(1:N/5) ) ;
xtitle(“signal blu bande supérieure”,“temps (s)”, “amplitude (V)”);
//
xsetech([0,3/4,1,1/4]) ; plot2d(t(1:N/5) ,s_dif(1:N/5) ) ;
xtitle(“différence des précédents”,“temps (s)”, “amplitude (V)”);
Modulations angulaires
**__Visualisation temporelle d’un signal modulé en phase__**
clear
// définition des constantes
N=400 ; fe=1e6 ; fp=50e3 ; finf=5e3; ind=5;
//
// description des vecteurs temps et signal
t=(0:N-1)/fe ;
//
//
sinf=1*sin(2*%pi*finf*t);
smod= 1*sin(2*%pi*fp*t + ind*sinf);
//
// initialisation de l'affichage
xbasc() ; xset( “font size”, 4);
//
// affichage du signal informatif
xsetech([0,0,1,1/2]) ; plot2d(t,sinf) ; xtitle(“signal informatif”,“temps (s)”,“amplitude (V)”);
//
// affichage du signal modulé
xsetech([0,1/2,1,1/2]) ; plot2d(t, smod) ;
xtitle(“signal modulé en phase”,“temps (s)”, “amplitude”);
clear
// constante de suréchantillage
N=16;
// vecteur de symboles
s_symb=[1 0 1 1 0 0 1 0];
// vecteur temps
t=(0 : 8*N-1);
// signal modulant
s_nrz=matrix(ones(N,1)*s_symb, 1, 8*N);
xbasc(); plot2d(t,s_nrz, rect=[0,-0.5,8*N,1.5]);
**__Signal modulé en fréquence__**
clear
//
// définition des constantes et du temps
N=1000 ; fe=1e6 ;f0=1e3; t=1/fe*(0:N-1);
//
// définition du vecteur fréquence
fp=f0*ones(1,N);
//
// définition de la phase (vecteur colonne)
for i=1:N,
theta(i)=2*%pi/fe*sum(fp(1:i));
end;
//
// calcul du sinus par les deux méthodes
//
s1=sin(theta');
s2=sin(2*%pi*f0*t);
//
// affichage
xbasc();
xsetech([0,0,1,1/2]) ; plot2d(t,s1) ; xtitle(“méthode de la phase”,“temps (s)”,“amplitude (V)”);
xsetech([0,1/2,1,1/2]) ; plot2d(t,s2) ; xtitle(“méthode classique”,“temps (s)”,“amplitude (V)”);
clear
//
// augmentation de la taille de la mémoire allouée à Scilab
stacksize(1.5e6);
//
// définition des constantes et du temps
N=1000 ; fe=1e6 ;f0=1e3; t=1/fe*(0:N-1);
//
// définition du vecteur fréquence
fp=f0*ones(1,N);
//
// définition du vecteur phase
theta=2*%pi/fe*fp*triu(ones(N, N));
//
// calcul du sinus par les deux méthodes
s1=sin(theta);
s2=sin(2*%pi*f0*t);
//
// affichage
xbasc();
xsetech([0,0,1,1/2]) ; plot2d(t,s1) ; xtitle(“méthode de la phase”,“temps (s)”,“amplitude (V)”);
xsetech([0,1/2,1,1/2]) ; plot2d(t,s2) ; xtitle(“méthode classique”,“temps (s)”,“amplitude (V)”);
clear; stacksize(1.5e6);
//
// définition des constantes et du temps
N=500 ; fe=1e6 ; fp=50e3 ; Finf=4e3; Uinf=2; t=1/fe*(0:N-1); ind=7;
//
// définition du signal informatif
sinf=Uinf*sin(2*%pi*Finf*t);
//
// définition du vecteur fréquence
f_inst=fp*ones(1,N)+ ind* Finf /Uinf*sinf;
//
// définition de la phase (vecteur colonne)
theta=2*%pi/fe*f_inst*triu(ones(N, N));
//
// calcul du signal modulé
smod=sin(theta);
//
// affichage
xbasc();xset(“font size”,4);
xsetech([0,0,1,1/2]) ; plot2d(t,sinf) ; xtitle(“signal modulant”,“temps (s)”,“amplitude (V)”);
xsetech([0,1/2,1,1/2]) ; plot2d(t,smod) ; xtitle(“signal modulé”,“temps (s)”,“amplitude (V)”);
**__Observation spectrale__**
clear; stacksize(1.5e6);
// définition des constantes et du temps
N=1000 ; fe=1e6 ; fp=50e3 ; Finf=2e3; Uinf=2; t=1/fe*(0:N-1); ind=5;
//
// définition du signal informatif et de la fréquence instantanée
sinf=Uinf*sin(2*%pi*Finf*t);f_inst=fp*ones(1,N)+ ind* Finf /Uinf*sinf;
//
// définition de la phase et calcul du signal modulé
theta=2*%pi/fe*f_inst*triu(ones(N, N));
smod=sin(theta');
//
// définition de l'échelle de fréquence et calcul du spectre
f=fe/N*(0:N-1);
S=fft(smod,-1);
//
// affichage
xbasc(); xset(“font size”,4);
xsetech([0,0,1,1/2]); plot2d3(f,1/N*abs(S));
xtitle(“spectre du signal modulé”,“fréquence (Hz)”,“amplitude (Vs)”);
xsetech([0,1/2,1,1/2]); plot2d3(f(35:65),1/N*abs(S(35:65)));
xtitle(“spectre zoomé”,“fréquence (Hz)”,“amplitude (Vs)”);
**__Annexe 4 : la programmation__**
// “tfcbi” calcul la fft, divise par le nombre d'échantillons,
// met sous forme -fe/2, fe/2, et génère le vecteur fréquence pour l'affichage
//
function[S, f]=tfcbi(s,fe)
N=length(s)
S1=1/N*fft(s,-1);
S=[S1(N/2+1:N), S1(1:N/2)];
f=fe/N*[-N/2 : N/2-1];
//
// “repet” répète un vecteur N fois (pour créer un signal périodique)
//
function[S]=repet(s,N)
S=matrix(s' *ones(1, N), 1 , N*length(s));
//
clear;
//
// chargement du fichier texte des fonctions
getf('mes_fonctions.txt');
//
// définition des constantes
a=10 ; N=1000 ;Te=1e-3 ;
//
// description du signal
//
s1=[ones(1, a), zeros(1, N-a)] ;
//
// appel de la fonction
[S1, f]=tfcbi(s1,1/Te);
//
// calcul de la transformée avec seulement la fonction fft (pour comparaison)
S=fft(s1,-1);
// affichage
xbasc(); xset(“font size”, 4);
xsetech([0,0,1,1/2]); plot2d(abs(S)); xtitle(“courbe générée par fft”,“échantillons”,“amplitude”);
xsetech([0,1/2,1,1/2]); plot2d(f,abs(S1));
xtitle(“transformée de Fourier”,“ fréquence (Hz)”,“amplitude (Vs)”);
**__Annexe 5 : polynômes et fonctions de transfert__**
| clear; // // définition des paramètres w0=10; z=.3; // // définition des pôles et zéro pole_d=-w0*(z+%i*sqrt(1-z | 2)); zero_n=0; // // définition du numérateur et dénominateur n=2*z/w0*poly([zero_n],'p'); d=1/w0 | 2*poly([ pole_d, pole_d' ],'p'); // // fonction de transfert h=n/d h = .06p —————– 1 + .06p + .01p2 |
|---|---|---|
| // // définition de la fonction de transfert h=n/d h = .06p —————– 1 + .06p + .01p2 |
–>roots(d)
ans =
! - 3. + 9.539392i !
! - 3. - 9.539392i !
//
// définition de la plage de variation
f=[.01:.01:100];
//
// calcul des points
h=freq(n,d,2*%i*%pi*f);
//
// affichage
xbasc(); xset(“font size”,4);
xsetech([0,0,1,1/2]); plot2d(f,20*log10(abs(h)), logflag=[“ln”]);
xtitle(“diagramme de Bode”,“fréquence (Hz)”,“ gain (dB)”);
xsetech([0,1/2,1,1/2]); plot2d(f,phasemag(h), logflag=[“ln”]);
xtitle(“”,“ fréquence (Hz)”,“phase en °”);
| clear; w0=10; z=.3; p=poly(0,'p'); h=(2*z*p/w0)/(1+2*z*p/w0+(p/w0) | 2) h = .06p —————– 1 + .06p + .01p2 |
|---|
clear;
w0=10; z=.3;
p=poly([ 0, 1],'p',’c’);
h=(2*z*p/w0)/(1+2*z*p/w0+(p/w0)^2)
h =
.06p
1 + .06p + .01p2
clear;
//
// définition de la variable p
p=poly(0,'p');
//
// définition de la fonction de transfert
h1=(1+0.1119*p^2)/(1+1.9234*p+1.808*p^2) ;
h2=(1+0.5788*p^2)/(1+0.4164*p+1.062*p^2) ;
h3=(1+0.8003*p^2)/(1+0.071*p+0.948*p^2) ;
h=h1*h2*h3 ;
//
// définition d’un système linéaire continu
h1=syslin('c',h);
//
// affichage de h1 pour des fréquences allant de 0,01 à 10 Hz par pas de 0,01
xbasc();xset(“font size”,4); bode(h1,.01,10,.01);
**__Annexe 6 : fonctions de Bessel__**
clear ; stacksize(1.5e6);
//
// définition du nombre points, de l'indice de modulation maximal, du nombre d'harmonique
N=500; Mmax=12; harmo=8;
//
// définition de l'angle d'intégration
theta=%pi/N*[0:N-1];
//
// définition de la plage de variation de l'indice de modulation
m=Mmax/N*[0:N-1];
//
// définition d'une matrice contenant les différents indices
J=zeros(harmo+1,N);
//
// initialisation de l'affichage
xbasc(); xset(“font size”,5);
xtitle (“fonction de Bessel”,“indice modulation”, “amplitude relative”);
//
// calcul et affichage des fonctions pour 8 fréquences harmoniques et la porteuse
for n=0:harmo,
J(n+1,:)= 1/N*(sum (cos(m'* sin(theta) - n*ones(N,1)*theta),'c'))';
plot2d(m,J(n+1,:), rect=[ 0, -0.4, 12, 1]);
end;
Traitement du signal avec Scilab : convolution et approche du filtrage numérique
clear ;
//
// définition de la fréquence d'échantillonnage, du nombre de points de calcul,
// du nombre d'échantillons du filtre, des vecteurs temps et échantillons
Fe=2 ; Te=1/Fe; N=300 ; Ne=50;
t=1/Fe*[0:N-1];n=[0:Ne-1];
//
// définition de la réponse impulsionnelle du filtre
f0=.02 ; w0=2*%pi*f0 ;
h=1/Fe*w0*exp(-w0*Te*n);
//
// définition des entrées
e1=ones(1,N);
e2=10*sin(w0*t);
//
// calcul des sorties et ajustement du nombre de points
s1=convol(e1,h); s1=s1(1:N);
s2=convol(e2,h); s2=s2(1:N);
//
// affichage sur la fenêtre 0
xset(“window”,0);xbasc(0); xset(“font size”,4);
xsetech([0 ,1/3,1,1/3]) ;plot2d(t,e1) ; ;plot2d(t,s1) ;
xtitle(“réponse à un échelon”);
xsetech([0 ,2/3,1,1/3]) ; ;plot2d(t,e2) ; ;plot2d(t,s2) ;
xtitle(“réponse sinusoïdale à la fréquence de coupure”);
xsetech([0 ,0,1,1/3]) ;plot2d3(n,h, rect=[-1,0,Ne,w0/Fe]);
xtitle(“réponse impulsionnelle”);
**__Approche de la transformation en z__**
**__Tracé de la réponse fréquentielle__**
// définition d’un polynôme en z à partir des coefficients de h
H=poly(h,'z','c');
//
// définition de la variation de la pulsation, qui n’a de signification que pour un intervalle 0 0,5
w=(0:.01:.5);
//
// calcul des valeurs particulières de H
ft=freq(H,1,exp(%i*w));
// affichage sur la fenêtre 1
xset(“window”,1);xbasc(1); xset(“font size”,4);
xsetech([0 ,0,1,1/2]) ;plot2d(w,abs(ft));
xtitle(“réponse fréquentielle sur échelles linéaires, normalisée en pulsation”);
xsetech([0 ,1/2,1,1/2]) ;
plot2d(w(2:length(w)),20*log10(abs(ft(2:length(w)))), logflag=[“ln”]) ;
xtitle(“réponse fréquentielle sur échelles semi-log, normalisées en pulsation”);
**__Annexe 1 : programmes des illustrations__**
clear ;
//
// définition des constantes, nombre de points, nombre d’impulsions élémentaires
// durée d’une impulsion, nombre d’échantillon sur une réponse élémentaire
N=200 ; Nrep=30; inc=4; Ne=80 ;
f0=.01 ; w0=2*%pi*f0 ;
n=[0:Ne-1];
//
// impulsion et réponse impulsionnelle
e0=ones(1,inc) ;
h0=w0*exp(-w0*n);
//
// définition des matrices où seront placé les impulsions et réponse
h=zeros(Nrep, N);
e=zeros(Nrep,N);
//
// initialisation de l’affichage pour la fenêtre 0
xset(“window”,0); xbasc(0) ; xset(“font size”, 4);
//
// boucle de calcul et d’affichage des impulsions et réponses
for i=0:inc:inc*Nrep-1,
h(i+1, : )=[zeros(1,i), h0 ,zeros(1,N-Ne-i)];
e(i+1, : )=[zeros(1,i), e0 ,zeros(1,N-inc-i)];
xsetech([0 ,0,1,1/3]) ;plot2d21) ;
end
//
// calcul de la sortie
s=sum(h,’r’);
//
// affichage de la sortie et des commentaires
xsetech([0 ,2/3,1,1/3]) ;plot2d(s ) ;
xtitle (“somme des réponses élémentaires”);
xsetech([0 ,0,1,1/3]) ; xtitle (“entrée décomposée en impulsions”);
xsetech([0 ,1/3,1,1/3]) ; xtitle (“réponse à chaque impulsion”);
//
// affichage des impulsions et réponses décalées sur la fenêtre 1
xset(“window”,1); xbasc(1) ; xset(“font size”, 4);
xsetech([0 ,0,1/2,1/2]) ;plot2d22) ;xtitle (“réponse h(n)”);
xsetech([1/2 , 0,1/2,1/2]) ;plot2d23) ;xtitle (“réponse h(n - k)”);
Traitement du signal avec Scilab : Signaux échantillonnés
**__Etude générale__**
clear
//
// définition des constantes, nombre de points par périodes de dirac, largeur du dirac, nombre de diracs, nombre de points total, fréquence et temps de calcul,
// fréquence et période du peigne de dirac et fréquence du signal
Npp=40 ; a=1; Nd=80; N=Npp*Nd;
Fe=16e3; Te=1/Fe; Fc=Fe*Npp; Tc=1/Fc; F=1000;
//
// description du vecteur temps
t=Tc*(0:N-1);
//
// signal à échantillonner
ve=5*sin(2*%pi*F*t);
//
// une période du dirac de rapport cyclique a/Npp
d=[ones(1, a), zeros(1, (Npp-a))];
//
// peigne de dirac sur Nd périodes
imp=matrix(d' *ones(1, Nd ), 1 , N);
//
// signal échantillonné
se=ve .* imp;
//
// réponse impulsionnelle du bloqueur
b=ones(1, Npp);
//
// signal en sortie du bloqueur, réajusté sur N points
v1= convol(b,se);
v1=v1(1 : N);
//
// calcul du filtre passe bas, butteworth d'ordre 8
G=analpf(8,'butt',[0 0],2*%pi*(0.4*Fe));
f=Fc/N*[0 : N-1];
GAIN=freq(G(2),G(3),%i*2*%pi*f);
G=[GAIN(1:N/2), conj(GAIN(N/2+1:-1:2))];
hpb=fft(G, 1);
//
// filtrage
v2=convol(hpb, v1); v2=v2(1:N);
// affichage
xset(“window”,0); xbasc(0) ; xset(“font size”,4);
xsetech([0,0,1,1/4]) ; plot2d(t,ve); plot2d2(t,imp);
xtitle(“signal entré, peigne de Dirac”,”temps (s)”,”amplitude (V)”)
xsetech([0,1/4,1,1/4]) ; plot2d2(t,se);
xtitle(“signal échantillonné”,”temps (s)”,”amplitude (V)”);
xsetech([0,2/4,1,1/4]) ; plot2d2(t,v1);
xtitle(“signal bloqué”,”temps (s)”,”amplitude (V)”)
xsetech([0,3/4,1,1/4]) ; plot2d2(t,v2);
xtitle(“signal de sortie filtré”,”temps (s)”,”amplitude (V)”);
**__Etude du modulateur__**
clear ;
// rappel des paramètres, redéfinition de l’amplitude de imp
Npp=40 ; a=20; Nd=20; N=Npp*Nd;
Fe=16e3; Te=1/Fe; Fc=Fe*Npp; Tc=1/Fc;
// description du vecteur temps
t=Tc*(0:N-1);
//
// une période du dirac de rapport cyclique a/Npp
d=[ones(1, a), zeros(1, (Npp-a))];
//
// peigne de dirac d’amplitude 1/a sur Nd périodes
tau=a*Te/Npp ;
imp=1/tau*matrix(d' *ones(1, Nd ), 1 , N);
//
// analyse de Fourier
f=Fc/N*[0:N/2-1];
IMP=1/N*abs(fft(imp,-1));
IMP=IMP(1:N/2);
//
// génération du sinus cardinal
S=abs(Fe*(sin(tau*%pi*f)./4)) +(f==0))) ;
//
// affichage
xset(“window”,0); xbasc(0) ; xset(“font size”,6);
xsetech([0,0,1,1/2]) ; plot2d2(t,imp, rect=([0, 0, 14E-4, 1.2*Npp/(a*Te)]) );
xtitle(“signal en fonction du temps”,”temps (s)”,”amplitude”)
xsetech([0,1/2,1,1/2]) ; plot2d(f,IMP, rect=([0, 0, 32e4, 1.2*Fe]));
plot2d(f(1 :4 :length(f)),S(1 :4 :length(f)), style=-3, rect=([0, 0, 32e4, 1.2*Fe]));
xtitle(“signal en fonction de la fréquence”,”fréquence (Hz) ”,”amplitude”);
// signal à échantillonner
F=3000 ; ve=5*sin(2*%pi*F*t);
//
// signal échantillonné
se=ve .* imp;
VE=1/N*abs(fft(ve,-1));
VE=VE(1:N/2);
//
SE=1/N*abs(fft(se,-1));
SE=SE(1:N/2);
//
// affichage
xset(“window”,1); xbasc(1) ; xset(“font size”,4);
xsetech([0,0,1,1/3]) ; plot2d(f, VE );
xtitle(“signal en entrée”,”fréquence”,”amplitude”)
xsetech([0,1/3,1,1/3]) ; plot2d(f, IMP );
xtitle(“impulsions”,”fréquence”,”amplitude”)
xsetech([0,2/3,1,1/3]) ; plot2d(f, SE);
xtitle(“signal échantillonné”,”fréquence”,”amplitude”);
**__Etude du démodulateur__**
clear ;
// rappel des paramètres
Npp=40 ; a=1; Nd=160; N=Npp*Nd; Fe=16e3; Te=1/Fe; Fc=Fe*Npp; Tc=1/Fc;
//
// description du vecteur temps et fréquence avec un « zoom » pour la fréquence
t=Tc*(0:N-1);
z=10 ; f=Fc/N*[0:N/z-1];
//
// peigne de dirac d’amplitude 1/a sur Nd périodes
tau=a*Te/Npp ;
d=[1*ones(1, a), zeros(1, (Npp-a))];imp=1/tau*matrix(d' *ones(1, Nd ), 1 , N);
//
// signal à échantillonner et échantillonnage
F1=2000 ;F2=4000 ; ve=5*(sin(2*%pi*F1*t)+sin(2*%pi*F2*t));
se=ve .* imp;
//
// réponse impulsionnelle du bloqueur
b=ones(1, Npp);b2=[b, zeros(1,N-Npp)];
//
// analyse de Fourier
IMP=1/N*abs(fft(imp,-1));IMP=IMP(1:N/z);
VE=1/N*abs(fft(ve,-1)); VE=VE(1:N/z);
SE=1/N*abs(fft(se,-1)); SE=SE(1:N/z);
B=1/N*abs(fft(b2,-1)); B=0.01*B(1 :N/z) ;
V1=SE.*B ;
//
// affichage
xset(“window”,0); xbasc(0) ; xset(“font size”,4);
xsetech([0,0,1,1/3]) ; plot2d(f, SE, rect=[0, 0, 7e4, 1.2*Fe*2.5] );
xtitle(“spectre du signal en entrée”,”fréquence”,”amplitude”)
xsetech([0,1/3,1,1/3]) ; plot2d(f, B );
xtitle(“réponse fréquentielle du bloqueur”,”fréquence”,”amplitude”)
xsetech([0,2/3,1,1/3]) ; plot2d(f, V1);
xtitle(“spectre du signal en sortie du bloqueur”,”fréquence”,”amplitude”);
// calcul du filtre et filtrage
G=analpf(8,'butt',[0 0],2*%pi*(0.4*Fe));
GAIN=freq(G(2),G(3),%i*2*%pi*f);
V2=GAIN.*V1 ;
//
// affichage
xset(“window”,1); xbasc(1) ; xset(“font size”,4);
xsetech([0,0,1,1/6]) ; plot2d(f, VE, rect=[0, 0, 4e4, 5] );
xtitle(“spectre du signal en entrée”,”fréquence”,”amplitude”)
xsetech([0,1/5,1,1/6]) ; plot2d(f, SE, rect=[0, 0, 4e4, 1.2*Fe*2.5] );
xtitle(“signal modulé”,”fréquence”,”amplitude”)
xsetech([0,2/5,1,1/6]) ; plot2d(f, abs(V1) , rect=[0, 0, 4e4, 5]);
xtitle(“spectre du signal en sortie du bloqueur”,”fréquence”,”amplitude”);
xsetech([0,3/5,1,1/6]) ; plot2d(f, abs(GAIN ) , rect=[0, 0, 4e4, 1.1]);
xtitle(“réponse fréquentielle du filtre”,”fréquence”,”amplitude”)
xsetech([0,4/5,1,1/6]) ; plot2d(f, abs(V2), rect=[0, 0, 4e4, 5] );
xtitle(“spectre du signal restitué”,”fréquence”,”amplitude”)
**__Annexe 1 : le sinus cardinal__**
clear ;
Delta=0.01 ; N=5*%pi ;
x=(0:Delta :N-1);
//
// génération du sinus cardinal
S=sin(x)./5)) +(x==0) ;
//
// affichage
xset(“window”,0); xbasc(0) ; xset(“font size”,6);
xsetech([0,0,1,1/2]) ; plot2d(x,S, axesflag=[5]);
xtitle(“sinc(x)=sin(x)/x”,””,””)
xsetech([0,1/2,1,1/2]) ; plot2d(x, S^2, axesflag=[5]);
xtitle(“sinus cardinal au carré”,””,””);
Traitement du signal avec Scilab : la transformée de Fourier discrète
**__Efficacité de l’algorithme de la fonction « fft »__**
clear ;
//
// nombre d’échantillons, 4 puis 512, puis1024 pour comparer les résultats
N=4 ;
//
// vecteur temps et signal
t=(0 :N-1) ; s=sin(t) ;
//
// début de la temporisation
timer() ;
//
// notre algorithme
Sf1=sum(s’ * ones(1,N) .* exp( -2*%pi * %i * t’ *t /N),’r’) ;
timer()
// lire la valeur retournée
//
// nouvelle temporisation pour l’algorithme fft cette fois
timer() ;
Sf2=fft(s,-1) ;
timer()
// lire la valeur et comparer à la précédente
**__Vérification expérimentale__**
clear
//
// paramètres de l’acquisition
// nombre de points, fréquence du signal, amplitude, fréquence d’échantillonnage
N=32; Fx=26e3 ; A=1 ; Fe=64e3 ;
//
// paramètre de l’échelle analogique : nombre de points
Na=1024 ;
//
// définition du vecteur temps et du vecteur fréquence normalisé à N
ta=N/(Na*Fe)*(0 :Na-1) ; ; k=N/Na*(0 : Na-1) ;
//
// définition du signal temporel
sa= A*exp(2*%i*%pi*Fx*ta) ;
//
// échelle des temps et fréquence numérique
tn=1/Fe*(0:N-1) ; fn=(0 :N-1) ;
//
// échantillonnage du signal temporel, calcul de la fft
sn=A*exp(2*%i*%pi*Fx*tn) ;
Sf=abs(fft(sn,-1));
//
// expression théorique de la DFT d’une exponentielle
x=k/N-Fx/Fe;
S= A*sin(N*x*%pi) . / (sin(x*%pi)+(x==0)) + (x==0)*N;
//
// précision sur l’échelle des fréquences des points affichables
pts=zeros(1, N) ;
//
// affichage
xset(“window”, 0) ; xbasc(0) ; xset (“font size”, 4) ;
xsetech([ 0, 0, 1, 1/2]);
plot2d(ta, real(sa), rect=[0, -1.1*A, N/Fe, 1.1*A], style=0) ;
plot2d3(tn,real(sn)) ;
xtitle(“signal analogique et échantillons temporels” ,“temps(s)”, “amplitude (V)”) ;
//
xsetech([ 0, 1/2, 1, 1/2]) ;
plot2d(k, abs(S), style=0, rect=[0, 0, N, 1.1*A*N]) ;
plot2d3(fn, Sf) ;
plot2d(fn, pts, style=-3) ;
xtitle(“transformée de Fourier théorique et calculée”,“rang”, “amplitude(V.s)”) ;
clear
// paramètre de l’acquisition
// nombre de points, fréquence du signal, amplitude
N=128; Fx=26e3 ; A=1 ; Fx2=27e3 ;
//
// paramètre de l’échelle analogique
// nombre de points, fréquence d’échantillonnage
Na=1024 ; Fe=64e3 ;
//
// définition du vecteur temps et du vecteur fréquence normalisé à N
ta=N/(Na*Fe)*(0 :Na-1) ; k=N/Na*(0 : Na-1) ;
//
// définition du signal temporel
sa= A*(exp(2*%i*%pi*Fx*ta)+ exp(2*%i*%pi*Fx2*ta)) ;
//
// échantillonnage du signal temporel, calcul de la fft
tn=1/Fe*(0:N-1) ; fn=(0 :N-1) ;
sn=A*(exp(2*%i*%pi*Fx*tn)+ exp(2*%i*%pi*Fx2*tn)) ;
Sf=abs(fft(sn,-1));
//
// expression théorique de la DFT d’une exponentielle
x=k/N-Fx/Fe;
x2=k/N-Fx2/Fe;
S= A*sin(N*x*%pi). / (sin(x*%pi)+(x==0)) + (x==0)*N;
S2= A*sin(N*x2*%pi). / (sin(x2*%pi)+(x2==0)) + (x2==0)*N;
// précision sur l’échelle des fréquences des points affichables
pts=zeros(1, N) ;
//
// affichage
xset(“window”, 0) ; xbasc(0) ; xset (“font size”, 4) ;
xsetech([ 0, 1/2, 1, 1/2]) ;
plot2d(k, abs(S), style=0, rect=[0, 0, N, 1.1*A*N]) ;
plot2d(k, abs(S2), style=0) ;
plot2d3(fn, Sf) ;
plot2d(fn, pts, style=-3) ;
xtitle(“transformée de Fourier calculée”,“rang”, “amplitude(V.s)”) ;
xsetech([ 0, 0, 1, 1/2]);
plot2d(ta, real(sa), style=0) ;
plot2d3(tn,real(sn)) ;
xtitle(“signal analogique et échantillons temporels” ,“temps(s)”, “amplitude (V)”) ;
//
xset(“window”, 1) ; xbasc(1) ; xset (“font size”, 4) ;
plot2d(k(0.33*Na :0.5*Na), abs(S(0.33*Na :0.5*Na)), rect=[0.33*Na/N, 0, 0.5*Na/N, 1.1*A*N] , style=0) ;
plot2d(k(0.33*Na :0.5*Na), abs(S2(0.33*Na :0.5*Na)) ,style=0) ;
plot2d3(fn(0.33*N :0.5*N), Sf(0.33*N :0.5*N)) ;
plot2d(fn(0.33*N :0.5*N), pts(0.33*N :0.5*N), style=-3) ;
xtitle(“transformée de Fourier calculée”,“rang”, “amplitude(V.s)”) ;
**__Annexe : programme des illustration__**
clear;
//
// définition des constantes, rapport cyclique, nombre de points affichés
// période et fréquence d’échantillonnage, amplitude du signal
// largeur de l’affichage (en nombre paire de périodes), constantes d’affichage
// vecteur temps
a=20 ; N=100 ;Te=10e-3 ; Fe=1/Te; Amp=5 ;
P=2; At=N*P*Te/2; Af=P*Fe/2;
t=Te*(-N/2: N/2-1);
//
// description du signal porte
s=Amp*[zeros(1,(N-a)/2), ones(1, a), zeros(1, (N-a)/2)] ;
//
// transformé de Fourier sous forme bi-lattérale et vecteur fréquence
S1=1/N*fft(s,-1);
S=[S1(N/2+1:N), S1(1:N/2)];
f=Fe/N*[-N/2 : N/2-1];
//
// affichage du signal et de sa transformée de Fourier
xset(“window”, 0) ; xbasc(0) ; xset (“font size”, 4) ;
xsetech([0,0,1,1/2]);
plot2d2(t, s, rect=[-At, 0, At, 1.2*Amp]);
xtitle(“signal temporel”,“temps (s)”,“amplitude (V)”);
xsetech([0, 1/2,1,1/2]);
plot2d(f,abs(S), rect=[-Af, 0,Af, 1.2*a *Te*Amp]);
xtitle(“transformée de Fourier”,“ fréquence (Hz)”,“amplitude (Vs)”);
//
// périodisation de la transformée de Fourier
// définition d’une nouvelle échelle de fréquence
S1=[S(N/2+1:N), S(1:N/2)];
S2=Fe*matrix(S1' *ones(1, P), 1 , N*P);
f2=Fe/N*[-P*N/2 : P*N/2-1];
//
//
//
// signal échantillonné, définition d’une nouvelle échelle de temps
se= Amp*[zeros(1,P*N/4+(N-a)/2), ones(1, a), zeros(1,P*N/4+(N-a)/2)] ;
t2=Te*(-P*N/2 : P*N/2-1);
//
// affichage du signal échantillonné et de sa transformée
xset(“window”, 1) ; xbasc(1) ; xset (“font size”, 4) ;
xsetech([0,0,1,1/2]);
plot2d3(t, s, rect=[-At, 0, At, 1.2*Amp]);
plot2d (t2, se, style=-9);
xtitle(“signal temporel”,“temps (s)”,“amplitude (V)”);
xsetech([0, 1/2,1,1/2]);
plot2d(f2,abs(S2), rect=[-Af, 0, Af, 1.2*a*Te*Fe*Amp]);
xtitle(“transformée de Fourier”,“ fréquence (Hz)”,“amplitude (Vs)”);
//
// effet de la troncature, le spectre d’origine est multiplié par N
S2=N*Te*S2;
//
// affichage du signal tronqué
xset(“window”, 2) ; xbasc(2) ; xset (“font size”, 4) ;
xsetech([0,0,1,1/2]);
plot2d3(t, s, rect=[-At, 0, At, 1.2*Amp]);
plot2d (t, s, style=-9);
xtitle(“signal temporel”,“temps (s)”,“amplitude (V)”);
xsetech([0, 1/2,1,1/2]);
plot2d(f2,abs(S2), rect=[-Af, 0, Af, 1.2*a*Te*N*Amp]);
xtitle(“transformée de Fourier”,“ fréquence (Hz)”,“amplitude (Vs)”);
//
// signal temporel périodisé
s1=[s(N/2+1:N), s(1:N/2)];
s2=N*Te*matrix(s1'*ones(1,P), 1, N*P);
//
// affichage du signal temporel périodique
xset(“window”, 3) ; xbasc(3) ; xset (“font size”, 4) ;
xsetech([0,0,1,1/2]);
plot2d3(t2, N/Fe*s2, rect=[-At, 0, At, 1.2*N/Fe*Amp]);
xtitle(“signal temporel”,“temps (s)”,“amplitude (V)”);
xsetech([0, 1/2,1,1/2]);
plot2d3(f2,abs(S2), rect=[-Af, 0, Af, 1.2*a*Te*N*Amp]);
xtitle(“transformée de Fourier”,“ fréquence (Hz)”,“amplitude (Vs)”);
Traitement du signal avec Scilab : puissance et densité spectrale de puissance
**__Application à un signal NRZ aléatoire__**
clear
//
// constantes : coefficient de sur-échantillonnage, nombre de symboles,
// nombre de points, période et fréquence d'échantillonnage, amplitude
cse=8; Nb=128; N=Nb*cse ; Te=1e-3; Fe=1/Te; A=5 ;
//
// vecteur temps de l'information et du signal
//
ind=(0:Nb-1);
t_nrz=Te*(0:N-1);
f=Fe/N*(0:N/2-1);
//
// signal aléatoire à + ou –1, répartition « normale »
s=sign(rand(1,Nb,'n'));
//
// sur-échantillonnage de s
s_nrz=A*(matrix(ones(cse,1)*s,1,N));
//
// calcul de la DSP
DSP=(1/N*abs(fft(s_nrz,-1)))^2;
DSPdB=10*log10( DSP(1:N/2) +%eps);
//
//
// affichage
xset(“window”, 0) ; xbasc(0) ; xset (“font size”, 4) ;
//
xsetech([ 0, 0, 1, 1/3]) ; plot2d3(ind, s, rect=[0, -1.2, Nb, 1.2]) ;
xtitle(“niveaux logiques émis” ,“rang”, “niveau”) ;
//
xsetech([ 0, 1/3, 1, 1/3]) ; plot2d2(t_nrz, s_nrz, rect=[0, -1.2*A, N*Te, 1.2*A]) ;
xtitle(“signal NRZ”,“temps (s)”,“amplitude (V)”) ;
//
xsetech([ 0, 2/3, 1, 1/3]) ; plot2d(f, DSPdB) ;
xtitle(“densité spectrale de puissance”, “fréquence (Hz)”, “DSP (dB)”) ;
//
//
// calcul de la puissance moyenne et vérification de Parseval
PMOY_T=1/N*sum(s_nrz^2)
PMOY_F=sum(DSP)
Traitement du signal avec Scilab : signaux aléatoires et bruits
**__Application__**
clear
//
// définition des constantes
loi='n';N=4096; V_moy=2; p=.5; x=(-6:.1:6);
//
// génération du bruit
bruit=sqrt(p)*rand(1,N,loi)+V_moy;
//
// calcul de la DSP
dsp=1/N*(abs(fft(bruit,-1)))^2 ;
//
// courbe théorique de probabilité
p_g=1/(sqrt(p)*sqrt(2*%pi))*exp(-6))
xtitle(“densité spectrale de puissance”,“fréquence”,“ DSP”) ;
//
// calcul des puissances
puissance_f=1/N*sum(dsp)
puissance_t=1/N*sum(bruit^2)
**__Annexe : programme des illustrations__**
clear;
//
// définition des constantes et variable d’ordonnées
mx=10 ; pas=0.1 ;
x=(-mx : pas : mx);
m=0 ; sig=2.5 ;
//
// variable gaussienne
p_g=1/(sig*sqrt(2*%pi)*exp(-7);
//
// variable uniforme
n=6 ;
p_u=1/n*[zeros(1,(mx-n/2)), ones(1,n), zeros(1, (mx-n/2))] ;
//
// variable discrète
nb=2 ;
p_disc=1/nb*[zeros(1,9), 1, 0, 1, zeros(1, 8)] ;
//
// affichage
xset(“window”,0) ; xbasc(0); xset(“font size”,4);
xsetech([ 0, 1/3, 1, 1/3]) ; ; plot2d(x,p_g)
xtitle(“loi de probabilité gaussienne”,“ valeur probable” ,“probabilité”) ;
//
xsetech([ 0, 0, 1, 1/3]) ; plot2d28);
xset(“window”,0) ; xbasc(0); xset(“font size”,4);
plot2d(x,p) ;
xtitle(“ probabilité gaussienne, moyenne -2, écart type 2,5”,“ valeur probable” ,“probabilité” );
Traitement du signal avec Scilab : corrélation
Réalisation pratique en numérique
clear ;
//
// définition du temps et des signaux
t=(0 :9) ;x=5*ones(1,10) ;y=x ;
//
// affichage
xbasc() ; xset(«font size»,4) ;
xsetech([0,0,1,1/2]) ; plot2d3(t,x,rect=[-10,0,15,6]) ;xtitle(« signal x ») ;
xsetech([0,1/2,1,1/2]) ; plot2d3(t,y,rect=[-10,0,15,6]) ;xtitle(« signal y »)
clear ;
//
// définition du temps et des signaux
t=(-10 :9) ;x=[zeros(1,10),5*ones(1,10)] ;y=[zeros(1,2),5*ones(1,10),zeros(1,8)] ;
//
// affichage
xbasc() ; xset(«font size»,4) ;
xsetech([0,0,1,1/2]) ; plot2d3(t,x,rect=[-10,0,15,6]) ;xtitle(« signal x ») ;
xsetech([0,1/2,1,1/2]) ; plot2d3(t,y,rect=[-10,0,15,6]) ;xtitle(« signal y décalé de 8 échantillons »)
Exemples de solutions
**__Théorème de Wiener-Kintchine pour un signal aléatoire__**
clear
//
//génération d’un signal aléatoire de puissance 3 V2
bruit=sqrt(3)*rand(1,4096,'n');
//
// calcul de la DSP avec pspect
[dsp1]=10*log10(pspect(128,256,'tr',bruit));
//
// calcul de l’autocorrélation et de sa transformation de Fourier
[cov]=corr(bruit,256);
dsp2=10*log10(abs(fft(cov,-1)));
//
//initialisation de l’affichage, réglage pour une taille de fonte de 4, création d’une variable t
xbasc(); xset(“font size”,4); tau=(0:255);
//
// affichage
xsetech([0,0,.5,.5]);plot2d(bruit);
xtitle(“représentation temporelle du bruit”);
//
xsetech([0,.5,.5,.5]);plot2d(dsp1);
xtitle(“densité spectrale de puissance”);
//
xsetech([.5,0,.5,.5]);plot2d(tau,cov,rect=[-20,-0.2,256,3.5]);
xtitle(“autocorrélation”);
//
xsetech([.5,.5,.5,.5]);plot2d(dsp2);
xtitle(“transformée de Fourier de l autocorrélation”)
Périodicité d’un signal NRZ pseudo-aléatoire
clear
//
//génération des informations binaires aléatoire
s=sign(rand(1,64,'n'));
//
// génération du signal NRZ à partir des informations binaires par suréchantillonnage.
s_nrz=5*(matrix(ones(64,1)*s,1,4096));
//
// autocorrélation
[cov]=corr(s_nrz,4096);
//
//affichage des résultats
t=(0:4095);tau=t;
xbasc(); xset(“font size”,4);
xsetech([0,0,1,1/3]);plot2d(t,s_nrz,rect=[0,-5.5,4095,5.5]);xtitle(“représentation temporelle”);
xsetech([0,1/3,1,1/3]);plot2d(tau,cov);xtitle(“autocorrélation”);
xsetech([0,2/3,1,1/3]);plot2d(tau,cov,rect=[0,-10,500,30]);xtitle(“autocorrélation dilatée”)
clear
//
//signal de base NRZ de 64 bits suréchantillonné d’un rapport 8
s=sign(rand(1,64,'n')); s_nrz=5*(matrix(ones(8,1)*s,1,512));
//
// 8 périodes de signal 64 bits
s_per=matrix(s_nrz'*ones(1,8),4096);
//
// autocorrélation
[cov_per]=corr(s_per,4096);
//
//affichage des résultats
t=(0:511);t_per=(0:4095); tau=t_per;
xbasc(); xset(“font size”,4);
//
xsetech([0,0,1,1/3]);plot2d(t,s_nrz,rect=[0,-5.5,511,5.5]);
xtitle(“signal de base”);
//
xsetech([0,1/3,1,1/3]);plot2d(t_per,s_per,rect=[0,-5.5,4095,5.5]);
xtitle(“signal complet”);
//
xsetech([0,2/3,1,1/3]);plot2d(tau,cov_per,rect=[-10,-10,4096,30]);
xtitle(“autocorrélation”);
Détection d’un signal noyé dans du bruit
clear;
// définition du temps
t=1e-6*(1:4096);
//
// définition du rapport signal sur bruit et des signaux
SB=.1; b=sqrt9);;
xtitle(“signal bruité”);
//
xsetech([0,2/3,1,1/3]); plot2d(cov);
xtitle(“autocorrélation”)
clear;
//
// définition du temps
t=1e-6*(1:4096);
//
// définition du rapport signal sur bruit et des signaux
SB=.001;
b=sqrt10);;xtitle(“signal bruité”);
xsetech([0,2/4,1,.22]); plot2d(t(1 :512),s_recep(1 :512), rect=[0,-6,512e-6,6]);;xtitle(“signal local”);
xsetech([0,3/4,1,.22]); plot2d(cov);xtitle(“corrélation du signal bruité et du signal local”)
Traitement du signal avec Scilab : transmission numérique en bande de base
Cas d’un canal non bruité de bande passante réduite
Aspect fréquentiel
clear
//
// définition des constantes
// période d’un symbole, coefficient de sur échantillonnage, fréquence de coupure du filtre
Ts=1e-3; Sech=16 ; fc=800 ;
// période d’échantillonnage, nombre de bits, nombre de points
Te=Ts/Sech ; Nbit=256 ;Nmax=Nbit*Sech ;
//
// définition du vecteur temps
t=Te*(0: Nmax-1);
//
// signal aléatoire, puis sur échantillonnage pour obtenir un signal NRZ
se=sign(rand(1,Nbit,'n'));
se_sur=(matrix(ones(Sech,1)*se,1,Nmax)) ;
//
// élaboration d’un filtre analogique de Chebycheff, d’ordre 4, d’ondulation 0,1
// de fréquence de coupure 800 Hz
Ge=analpf(4,'cheb1',[.1 0],2*%pi*fc);
//
// préparation du tracé de la réponse fréquentielle du filtre
f=1/Te/Nmax*(1:Nmax);
gain=freq(Ge(2),Ge(3),%i*2*%pi*f);
//
// calcul des DSP du signal NRZ et signal NRZ filtré
DSPe=(1/Nmax*abs(fft(se_sur,-1)))^2;
DSPr=DSPe.*abs(gain).^2;
//
//
//
//
// troncature de grandeurs à afficher pour effet zoom
f_aff=f(.11/16*Nmax:Nmax/4);
gain_aff=abs(gain(.11/16*Nmax:Nmax/4));
DSPe_aff=DSPe(.11/16*Nmax:Nmax/4);
DSPr_aff=DSPr(.11/16*Nmax:Nmax/4);
//
//affichage fréquentiel
xset(“window”,0) ; xbasc(0); xset(“font size”,4); xset(“thickness”,2);
xsetech([0,0,1,.3]);
plot2d(f_aff, 10*log10(DSPe_aff+%eps),logflag=“ln”,rect=[1e2,-180,1e4,20],nax=[1,2,2,5]);
xtitle(“DSP du signal transmis”,“f(Hz)”,“DSP”);
xsetech([0,1/3,1,.3]);
plot2d(f_aff,20*log10(gain_aff),logflag=“ln”, rect=[1e2,-70,1e4,10],nax=[1,2,2,4]);
xtitle(“gain du filtre”,“f(Hz)”,“DSP”);
xsetech([0,2/3,1,.3]);
plot2d(f_aff, 10*log10(DSPr_aff+%eps),logflag=“ln”, rect=[1e2,-180,1e4,20],nax=[1,2,2,5]);
xtitle(“DSP du signal reçu”,“f(Hz)”,“DSP”);
Aspect temporel
// calcul de la fft et filtrage
Se_sur=1/Nmax*fft(se_sur,-1) ;
Sr=Se_sur .* gain ;
//
// calcul de la fft inverse en respectant les symétrie
Sr=[Sr(1 :Nmax/2), conj(Sr(Nmax/2+1 : -1 : 2))];
sr= Nmax*fft(Sr, 1) ;
//
// affichage
xset(“window”,1) ; xbasc(); xset(“font size”,4) ; xset(“thickness”,2) ;
//
xsetech([0,0,1,1/3]);plot2d311),se,rect=[0,-1.5, Nbit+1,1.5]);
xtitle(“information émise”,“rang”, “amplitude (V)”);
//
xsetech([0,1/3,1,1/3]);plot2d(t, se_sur, rect=[0, -1.5, Nmax*Te, 1.5]);
xtitle(“signal émis”,“temps (s)”, “amplitude (V)”);
//
//
//
xsetech([0,2/3,1,1/3]);plot2d(t,sr, rect=[0, -1.5, Nmax*Te, 1.5]) ;
xtitle(“signal reçu”,“temps (s)”, “amplitude (V)”);
// troncature des signaux
se2=se(1:8);
se_sur2=se_sur(1:8*Sech);
sr2=sr(1:8*Sech);
t2=t(1:8*Sech);
//
// affichage
xset(“window”,2) ; xbasc(2) ; xset(“font size”,4) ; xset(“thickness”,2);
//
xsetech([0,0,1,1/3]);plot2d312),se2,rect=[0,-1.5,8+1,1.5]);
xtitle(“information émise”,“rang”, “amplitude (V)”);
//
xsetech([0,1/3,1,1/3]);plot2d2(t2,se_sur2, rect=[0, -1.5, 8*Sech*Te, 1.5]);
xtitle(“signal émis”,“temps (s)”, “amplitude (V)”);
//
xsetech([0,2/3,1,1/3]);plot2d(t2,sr2, rect=[0, -1.5, 8*Sech*Te, 1.5]) ;
xtitle(“signal reçu”,“temps (s)”, “amplitude (V)”);
Diagramme de l’œil
// Diagramme de l’œil
//
// mise en forme de la matrice du signal
s_oeil=matrix(sr,2*Sech,Nbit/2);
//
// mise en forme de la matrice du temps
t_oeil=Te*(1:2*Sech)' *ones(1,Nbit/2);
//
// mise en place des styles, tous identiques
style=ones(1,Nbit/2);
//
//
// affichage
xset(“window”,3) ; xbasc(3) ; xset(“font size”,5) ; xset(“thickness”,1);
plot2d(t_oeil, s_oeil,style,) ;
xtitle(“diagramme de l œil”,“temps (s)”, “amplitude (V)”);
Suppression de l’interférence entre symbole
// Réponse impulsionnelle du filtre
clear
//
// définition des constantes, dont le facteur d'arrondissement et nombre d'échantillons du filtre
Ts=1e-3; Sech=16 ; Te=Ts/Sech; r=.6; Ne=128;
//
// définition du vecteur temps pour 128 échantillons
t=Te*(-Ne/2:Ne/2);
//
// calcul de la réponse impulsionnelle du filtre
h1=cos(%pi*r*t/Ts)./(1-(4*r^2*t.^2)/Ts^2+(abs(r*t)==Ts/2))+(abs(r*t)==Ts/2)*%pi/4;
h2=(sin(%pi*t/Ts))./(%pi*t/Ts+(t==0))+(t==0);
h=h1.*h2;
//
// affichage
xbasc(); xset(“font size”,5); xset(“thickness”,2);
plot2d3(t,h);
xtitle(“réponse impulsionnelle pour un filtre en cosinus surélevé; r=0,6”, “temps (s)”, “amplitude”);
// Réponse fréquentielle
clear
//
// définition des constantes
Ts=1e-3; Nbpt=4096; pas=1/Ts/Nbpt; r=.6;
//
// définition du vecteur fréquence
f=pas*(1:Nbpt);
//
// calcul par zones de la réponse
z1=Ts*ones(1,Nbpt*(1-r)/2 -1);
z2=Ts/2*(1-(sin(%pi*Ts*(f-1/2/Ts)/(r+(r==0)))));
z22=z2(Nbpt*(1-r)/2:Nbpt*(1+r)/2);
z3=zeros(1, Nbpt*(1-(1+r)/2)+1) ;
H=[z1,z22,z3];
//
// affichage
xbasc(); xset(“font size”,5); xset(“thickness”,2);
plot2d(f,H(1:Nbpt))
xtitle(“réponse fréquentielle du filtre en cosinus surélevé pour r=0,6 et Ts=1 ms”,“fréquence (Hz)”);
// Réponses impulsionnelles paramètrées par r
clear
//
//
//
// définition des constantes
Ts=1e-3; Sech=16 ; Te=Ts/Sech; r=(0:.5:1)'; Nbit=256 ; Ne=128 ;
//
// vecteur temps
t=Te*(-Ne/2: Ne/2);
//
// calcul des réponses
h1=cos(%pi*r*t/Ts)./(1-(4*r^2*t.^2)/Ts^2+(abs(r*t)==Ts/2))+(abs(r*t)==Ts/2)*%pi/4;
h2a=(sin(%pi*t/Ts))./(%pi*t/Ts+(t==0))+(t==0);
h2=ones(3,1)*h2a;
h=h1.*h2;
//
// affichage
h_aff=h';
style=(1:3);
t_aff=t'*ones(1,3);
xbasc(); xset(“font size”,4); xset(“thickness”,2);
plot2d(t_aff,h_aff,style);
xtitle(“réponse impulsionnelle des filtres en cosinus surélevé paramétrés en r”);
// Réponse impulsionnelle du filtre
clear ;
//
// constantes
Ts=1e-3; Sech=8 ; Nech=64; Te=Ts/Sech;r=.6 ; Nbit=256
//
// vecteur temps
t=Te*(1:4096) ;
//
// réponse impulsionnelle du filtre
Tf=Te*(-Nech/2:(Nech/2-1)) ;
h1=cos(%pi*r*Tf/Ts)./(1-(4*r^2*Tf.^2)/Ts^2+(abs(r*Tf)==Ts/2))+(abs(r*Tf)==Ts/2)*%pi/4;
h2=(sin(%pi*Tf/Ts))./(%pi*Tf/Ts+(Tf==0))+(Tf==0);
h=h1.*h2;
//
// affichage de la réponse
xset(“window”,0) ; xbasc(0); xset(“font size”,5); xset(“thickness”,2);
plot2d313),h);
xtitle(“réponse impulsionnelle du filtres pour r=0,6 sur 64 échantillons”,“rang”, “amplitude”);
// Signaux obtenus
//
// signaux binaire et NRZ
se=sign(rand(1,Nbit,'n'));
se_sur=(matrix([1; zeros(Sech-1,1)]*se,1,Sech*Nbit)) ;
//
// signal filtré reçu par le récepteur
sr=convol(h,se_sur);
//
// suppression des échantillons introduit par la convolution
sr_aff=sr(1:length(sr)-Nech+1);
//
// dilatation des échelles d’affichage
se2=se(1:16);
se_sur2=se_sur(1:16*Sech);
sr_aff2=sr_aff(1+Nech/2:16*Sech+Nech/2);
t2=t(1:16*Sech);
//
// affichage
xset(“window”,1) ; xbasc(1); xset(“font size”,4); xset(“thickness”,2);
//
xsetech([0,0,1,1/3]);plot2d314),se2,rect=[0,-1.5,16+1,1.5]);
xtitle(“information émise”,“rang”,“niveau”);
//
xsetech([0,1/3,1,1/3]);plot2d2(t2,se_sur2, rect=[0,-1.5,16e-3,1.5]);
xtitle(“signal émis”,“temps (s)”,“amplitude (V)”);
//
xsetech([0,2/3,1,1/3]);plot2d(t2,sr_aff2, rect=[0,-1.5,16e-3,1.5]) ;
xtitle(“signal reçu”,“temps (s)”, “amplitude (V)”);
// Génération du signal NRZ par convolution
//
// réponse impulsionnelle du filtre
h_symb=ones(1,Sech);
//
// convolution
s_symb=convol(se_sur2,h_symb);
//
// affichage
xset(“window”,2) ; xbasc(2) ; xset(“font size”,4); xset(“thickness”,2);
xsetech([0,0,1,1/2]);plot2d315),se2,rect=[0,-1.5,16+1,1.5]);
xtitle(“information émise”,“rang”,“niveau”);
//
xsetech([0,1/2,1,1/2]); plot2d2(t2,s_symb(1:length(s_symb)-8+1)) ;
xtitle(“signal codé”, “temps (s)”,“amplitude (V)”);
// Diagramme de l’œil
//
// suppression des effets de bord
sr_rec=sr(1+Nech/2:length(sr)-Nech/2+1-2*Sech);
//
// mise en forme de la matrice
s_oeil=matrix(sr_rec,2*Sech,Nbit/2-1);
//
// matrice du temps de du style
t_oeil=(1:2*Sech)'*ones(1,Nbit/2-1);
style=ones(1,Nbit/2-1);
//
// affichage
xset(“window”,3) ;xbasc(3); xset(“font size”,4); xset(“thickness”,2);
plot2d(t_oeil, s_oeil,style);xtitle(“diagramme de l oeil”);
Cas d’un canal bruité de bande passante infinie
// Réponse à une impulsion et à un symbole du filtre adapté
clear;
//
// constantes, débit, facteur de sur échantillonnage, période d’échantillonnage
// nombres de bits émis, nombre de points
// rapport signal sur bruit, constante du temps du filtre adapté
Ts=1e-3; Sech=16 ; Te=Ts/Sech;
Nbit=2; Nbpt=Nbit*Sech;
SB=100; t0=3*Ts/2;
//
// vecteurs temps
t=Te*(1:Nbpt);
//
// donnée 1 suivie de 0, puis signal NRZ correspondant
se=[1 0];
se_sur=(matrix(ones(Sech,1)*se,1,Nbpt)) ;
//
// réponse du filtre adapté
gr=[zeros(1,(-Ts+t0)/Te),ones(1,Ts/Te)];
//
// signal en sortie du filtre
sf=convol(gr,se_sur);
//
// adaptation du temps pour l’affichage complet du signal filtré
tsf=Te*(1:length(sf));
//
// affichage
xbasc(); xset(“font size”,4); xset(“thickness”,2);
//
xsetech([0,0,1,1/3]); plot2d(t,se_sur,rect=[0,0,2e-3,1.5]);
xtitle(“signal émis”,“ temps”, “amplitude”);
//
xsetech([0,1/3,1,1/3]);plot2d([1:length(gr)],gr,rect=[0,0,2*Sech,1.5]);
xtitle(“réponse impultionnelle du filtre”,“rang”, “amplitude”);
//
xsetech([0,2/3,1,1/3]); plot2d(tsf,sf);
xtitle(“réponse du filtre au signal émis”,“temps”, “amplitude”);
// Réponse d’une chaîne complète avec filtre adapté
//
clear;
//
// constantes
Ts=1e-3; Sech=8 ; Te=Ts/Sech; Nbit=16; Nbpt=Nbit*Sech; t=Te*(1:Nbpt);SB=100; t0=3*Ts/2;
//
// données et signal NRZ correspondant
se=sign(rand(1,Nbit,'n'))
se_sur=(matrix(ones(Sech,1)*se,1,Nbpt)) ;
//
// signal fictif permettant de repérer les instants idéaux d’échantillonnage et les valeurs attendues
// suréchantillonnage par des 0 du signal de donnée
se_surd=(matrix([1; zeros(Sech-1,1)]*se,1,Nbpt)) ;
//
// synchronisation avec les maxima et minima de la réponse et troncature pour l’affichage
ck=3*[zeros(1,t0/Te-1),se_surd]; ck=ck(1: Nbpt);
//
// génération du bruit
b=sqrt(1/SB)*rand(1,Nbpt,'n');
//
//
// ajout du bruit au signal utile
sr=se_sur+b;
//
// réponse du filtre
gr=[zeros(1,(-Ts+t0)/Te),ones(1,Ts/Te)];
//
// filtrage et troncature pour l’affichage
sf=convol(gr,sr); sf=sf(1: Nbpt);
//
// extraction du signe de la valeur au moment de l’échantillonnage, donc du niveau logique.
sortie=sign(sf.*abs(ck)) ;
//
// mise sous forme NRZ par convolution de la valeur et troncature pour affichage
s_basc=convol(ones(1,Sech), sortie); s_basc=s_basc(1:Nbpt);
//
// affichage
xbasc(); xset(“font size”,3); xset(“thickness”,2);
//
xsetech([0,0,1,1/4]);plot2d316),se,rect=[0,-1.5,Nbit+1,1.5]);
xtitle(“information émise”);
//
xsetech([0,1/4,1,1/4]);plot2d(t,sr); xtitle(“signal bruité reçu”);
//
xsetech([0,2/4,1,1/4]);plot2d(t,sf);xtitle(“signal en sortie du filtre et instants de décision”);
xsetech([0,2/4,1,1/4]);plot2d3(t,ck/2);
//
xsetech([0,3/4,1,1/4]);plot2d3(t,ck/8,rect=[0,-1.5,16e-3,1.5]);
//
xsetech([0,3/4,1,1/4]);plot2d(t,s_basc);xtitle(“signal en sortie du récepteur”);
Traitement du signal avec Scilab : analyse et synthèse des filtres numériques
Analyse de quelques filtres
**__Filtre passe bas__**
clear
//
// constante
a=.5;
//
// définition du numérateur
h1=poly([1 1],'z','c');
//
// définition du dénominateur
h2=poly([-a 1],'z','c');
//
// gamme de fréquence normalisée par rapport à la fréquence d’échantillonnage
f=(0:.01:1);
//
// calcul des différents points de la fonction de transfert
hf=freq(h1,h2,exp(2*%pi*%i*f));
//
// affichage
xbasc(); xset(“font size”,4);
xsetech([0,0,1,.5]); plot2d(f,abs(hf))
xtitle(“module du filtre en fonction de la fréquence normalisée”)
xsetech([0,.5,1,.5]); plot2d(f,atan(imag(hf),real(hf)));
xtitle(“phase du filtre en fonction de la fréquence normalisée”)
// gamme de fréquence normalisée par rapport à la fréquence d’échantillonnage
f=(0:.01:.5);
//
// calcul des différents points de la fonction de transfert
hf=freq(h1,h2,exp(2*%pi*%i*f));
//
// affichage
xbasc(); xset(“font size”,4);
xsetech([0,0,1,.5]); plot2d(f,abs(hf));xtitle(“module du filtre en fonction de la fréquence normalisée”)
xsetech([0,.5,1,.5]); plot2d(f,atan(imag(hf),real(hf)));
xtitle(“phase du filtre en fonction de la fréquence normalisée”)
clear
//
// constantes et vecteur temps
a=.5; TE=1e-3;
t=TE*(0:127);
//
// définition d'une entrée
x=5*sin(2*%pi*0.01/TE*t)+3*sin(2*%pi*0.4/TE*t);
//
// définition de la fonction en z
h=poly([1 1],'z','c')./poly([-a 1],'z','c');
//
// création d’un système linéaire
hz=syslin('d',h);
//
// filtrage
y=flts(x,hz);
// affichage
xbasc(); xset(“font size”,4);
xsetech([0,0,1,.5]); plot2d(t,x);xtitle(“signal incident”,“temps (s)”,“amplitude”)
xsetech([0,.5,1,.5]); plot2d(t,y);xtitle(“signal de sortie”,“temps (s)”,“amplitude”)
Synthèse par la méthode des fenêtres
clear ;
//
// définition des constantes
Te=1e-3 ; Nb=15; Fc=.25/Te;
//
// paramètres “temps”
t=Te*(-(Nb-1)/2 : (Nb-1)/2);
//
// réponse impulsionnelle
g=2*Fc*Te*(sin(2*%pi*t*Fc)./(2*%pi*t*Fc+(t==0))+(t==0));
//
// calcul de la réponse fréquentielle à partir le transformée en z
f=(0:.001:.5);
h=poly(g,'z','c');
H=freq(h,1,exp(2*%pi*%i*f));
//
// affichage
xbasc(); xset(“font size”,4);
xsetech([0,0,1,.3]); plot2d(t,g);
xtitle(“réponse impulsionnelle du filtre idéal”,“temps”,“amplitude”) ;
//
xsetech([0,1/3 ,1,.3]); plot2d317),g);
xtitle(“réponse discrétisée et causale”,“échantillon”,“amplitude”) ;
//
xsetech([0,2/3 ,1,.3]); plot2d(f,abs(H));
xtitle(“réponse fréquentielle”,“f normalisée”,“gain linéaire”) ;
**__Linéarité de la phase__**
clear ;
//
// définition des constantes
Te=1e-3 ; Nb=15; Fc=.25/Te;
//
// paramètres “temps”
t=Te*(-(Nb-1)/2 : (Nb-1)/2);
//
// réponse impulsionnelle
g=2*Fc*Te*(sin(2*%pi*t*Fc)./(2*%pi*t*Fc+(t==0))+(t==0));
//
// calcul de la réponse fréquentielle à partir le transformée en z
f=(-1:.001:1);
h=poly(g,'z','c');
H=freq(h,1,exp(2*%pi*%i*f));
//
// affichage
xbasc(); xset(“font size”,4);
xsetech([0,0,1,.5]); plot2d(f,abs(H));xtitle(“réponse fréquentielle”,“f normalisée”,“gain linéaire”)
xsetech([0,.5 ,1,.5]); plot2d(f,atan(imag(H), real(H)));
xtitle(“phase”,“f normalisée”,“phase en rd”)
clear ;
//
// définition des constantes
Te=1e-3 ; Nb=16; Fc=.25/Te;
//
// paramètres “temps”
t=Te*(-(Nb-1)/2 : (Nb-1)/2);
//
// réponse impulsionnelle
g=2*Fc*Te*(sin(2*%pi*t*Fc)./(2*%pi*t*Fc+(t==0))+(t==0));
//
// calcul de la réponse fréquentielle à partir le transformée en z
f=(-1:.001:1);
h=poly(g,'z','c');
H=freq(h,1,exp(2*%pi*%i*f));
//
// affichage
xbasc(); xset(“font size”,4);
xsetech([0,0,1,.5]); plot2d(f,abs(H));xtitle(“réponse fréquentielle”,“f normalisée”,“gain linéaire”)
xsetech([0,.5 ,1,.5]); plot2d(f,atan(imag(H), real(H)));
xtitle(“phase”,“f normalisée”,“phase en rd”)
**__Choix d’une fenêtre__**
clear ;
//
// fenêtre rectangulaire
[hrc,hrm,fr]=wfir('lp',33,[.25 .4],'re',[0 0]);
//
// fenêtre de Kaiser avec b=5,6
[hkc,hkm,fk]=wfir('lp',33,[.25 .4],'kr',[5.6 0]) ;
//
// fenêtre de Hamming avec a=0,54
[hhc,hhm,fh]=wfir('lp',33,[.25 .4],'hm',[0.54 0]) ;
//
// affichage
xbasc() ; xset (“font size”,4);
xsetech([0,0,1,.3]); plot2d(fr,hrm);xtitle(“fenêtre rectangulaire”);
xsetech([0,1/3,1,.3]); plot2d(fk,hkm);xtitle(“fenêtre de Kaiser avec 5,6”);
xsetech([0,2/3,1,.3]); plot2d(fh,hhm);xtitle(“fenêtre de Hamming avec 0,54”);
Méthode de l’échantillonnage en fréquence
clear
//
// gabarit analogique synthétisé du filtre
Ha=[0*ones(1,15) ones(1,10) 0*ones(1,39)];
//
// calcul de la réponse fréquentielle réelle
Hd=fsfirlin(Ha,1);
//
// paramètres d'affichage
fd=0.5/length(Hd)*(0:length(Hd)-1);
fa=0.5/length(Ha)*(0:length(Ha)-1);
//
// affichage
xbasc();xset(“font size”,4);
plot2d(fa,abs(Ha),style=-2);
plot2d(fd,abs(Hd));
xtitle(“réponse du filtre numérique et points d échantillonnage du gabarit”);
Synthèse des filtres récursifs à réponse impulsionnelle infinie
**__Méthode de l’invariance impulsionnelle__**
//
clear
//
// définition des constante
Te=1e-3; w0=.01*2*%pi/Te; A=1;
//
// variable temps
t=Te*(0:128);
//
// réponse impulsionnelle
h=w0*exp(-w0*t);
//
// description du filtre par la transformée en z, numérateur puis dénominateur
H1n=A*w0*poly([0 1],'z','c');H1d= poly([-exp(-w0*Te) 1],'z','c');
//
// description d’un second filtre en ajustant le gain basse fréquence
A=(1-exp(-Te*w0))/w0;
//
H2n=A*w0*poly([0 1],'z','c');H2d= poly([-exp(-w0*Te) 1],'z','c');
//
//
//
// calcul de la réponse fréquentielle des deux filtres précédents
f=(0:.001:.05);
hdf1=freq(H1n,H1d,exp(2*%pi*%i*f));
hdf2=freq(H2n,H2d,exp(2*%pi*%i*f));
//
// affichage
xbasc(); xset(“font size”,4);
xsetech([0,0,1,.3]);plot2d(t,h, style=-9);plot2d3(t,h);xtitle(“réponse impulsionnelle”);
xsetech([0,1/3,1,.3]);plot2d(f,hdf1);xtitle(“réponse fréquentielle pour A=1”);
xsetech([0,2/3,1,.3]);plot2d(f,hdf2);xtitle(“réponse fréquentielle corrigée”);
Transformation bilinéaire
clear
//
// définition des constantes
Te=1e-3; Fcd=200;
//
// filtre analogique équivalent
Ga=analpf(6,'cheb1',[.1,0],2*%pi*Fcd);
//
// filtre analogique pour le calcul
teta=Fcd*(2*Te)*%pi;
Fcad=(2/Te*tan(teta/2))/(2*%pi);
Gad=analpf(6,'cheb1',[.1,0],2*%pi*Fcad);
//
// définition de la variable z
z=poly(0,'z');
//
// transformation bilinéaire
Gd=horner(Gad,2/Te*(z-1)/(z+1));
//
// calcul des points du filtre analogique équivalent
fa=(100:10:300);
Gain_a=freq(Ga(2),Ga(3),%i*2*%pi*fa);
//
// calcul des points du filtre numérique
fd=(.1:.01:.3);
Gain_d=freq(Gd(2),Gd(3),exp(%i*2*%pi*fd));
//
// affichage
xbasc(); xset(“font size”,4);
//
xsetech([0,0,1,.5]);plot2d(fa,20*log10(abs(Gain_a)), rect=[100, -60 300, 10]);
xtitle(“réponse fréquentielle du filtre analogique équivalent”,“fréquence”, “gain en dB”);
//
xsetech([0,1/2,1,.5]);plot2d(fd,20*log10(abs(Gain_d)) ,rect=[0.1, -60 0.3, 10]);
xtitle(“réponse fréquentielle du filtre numérique”,“f normalisée”,“gain en dB”);