1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131
|
// -------------------------------------------------------------------------
function [vect_filt]=filt_pb_iir(vect_donnees,freq_coupure,freq_echant,ordre)
Z = vect_donnees;
vect_filt = Z ;
N = length ( Z ) ;
dim = size ( Z ) ;
moy_Z = mean ( Z ) ;
Z = Z - moy_Z ;
// Recherche d'un IIR stable
[nu_c_est, num, den, bStable] = filnum_RII_pb(ordre,freq_coupure/freq_echant,'butt',-0.1,-40,512)
while ( bStable == %F ) & ( ordre > 0 )
ordre = ordre - 1 ;
[nu_c_est, num, den, bStable] = filnum_RII_pb(ordre,freq_coupure/freq_echant,'butt',-0.1,-40,512)
end;
if ( bStable ) & ( N == max ( dim(1) , dim(2) ) ) then
// ----- Filtrage dans le domaine temporel -----
// Filtrage dans le sens directe
Z_filt_f = zeros ( N , 1 ) ;
for i=(ordre+1):N
Z_in = Z ( i:-1:i-ordre ) ;
Y_in = Z_filt_f ( i:-1:i-ordre ) ;
Z_filt_f ( i ) = num * Z_in - den * Y_in ;
end ;
// Filtrage dans le sens retrograde
Z_filt_b = zeros ( N , 1 ) ;
for i=N-ordre:-1:1
Z_in = Z_filt_f ( i:i+ordre ) ;
Y_in = Z_filt_b ( i:i+ordre ) ;
Z_filt_b ( i ) = num * Z_in - den * Y_in ;
end ;
vect_filt ( : ) = Z_filt_b ( : ) + moy_Z ;
end;
endfunction
// -------------------------------------------------------------------------
function [nu_c_est, num, den, bStable] = filnum_RII_pb(ordre,nu_c,meth,GdB_min_bp, GdB_max_bc,N)
// entrées :
// ordre : ordre des polynômes numérateur et dénominateur
// nu_c : fréquence de coupure
// meth = 'cheb1', 'cheb2','ellip','butt'
// GdB_min_bp : Gain en dB minimum dans la bande passante - utile uniquement pour 'cheb1' et 'ellip'
// GdB_max_bc : Gain en dB maximum dans la bande coupée - utile uniquement pour 'cheb2' et 'ellip'
// sorties :
// nu_c_est : fréquence de coupure estimée à partir des N points
// num : coefficients du polynôme numérateur (puissance décroissante)
// den : coefficients du polynôme dénominateur (puissance décroissante)
// bStable : booleen permettant de savoir si le filtre est stable au sens EBSB ou non
gain_min_bp = 1-10^(GdB_min_bp/20);
gain_max_bc = 10^(GdB_max_bc/20);
for r = 1 : N
if((r-1)/(2*N) < nu_c)
Gab_id(r) = 1;
else
Gab_id(r) = 0;
end
D_m3dB(r) = -3;
end
if meth == 'cheb1' then
H = iir(ordre,'lp',meth,[nu_c 0], [gain_min_bp 0]);
str1 = 'Gain d''une fonction de transfert d''un filtre passe-bas RII obtenue avec la méthode de Chebyshev 1';
str2 = 'Gain en dB d''une fonction de transfert d''un filtre passe-bas RII obtenue avec la méthode de Chebyshev 1';
elseif meth == 'cheb2' then
H = iir(ordre,'lp',meth,[nu_c 0], [0 gain_max_bc]);
str1 = 'Gain d''une fonction de transfert d''un filtre passe-bas RII obtenue avec la méthode de Chebyshev 2';
str2 = 'Gain en dB d''une fonction de transfert d''un filtre passe-bas RII obtenue avec la méthode de Chebyshev 2';
elseif meth == 'ellip' then
H = iir(ordre,'lp',meth,[nu_c 0], [gain_min_bp gain_max_bc]);
str1 = 'Gain d''une fonction de transfert d''un filtre passe-bas RII obtenue avec la méthode elliptique';
str2 = 'Gain en dB d''une fonction de transfert d''un filtre passe-bas RII obtenue avec la méthode elliptique';
else
H = iir(ordre,'lp','butt',[nu_c 0], [0 0]);
str1 = 'Gain d''une fonction de transfert d''un filtre passe-bas RII obtenue avec la méthode de Butterworth';
str2 = 'Gain en dB d''une fonction de transfert d''un filtre passe-bas RII obtenue avec la méthode de Butterworth';
end
[H_m,nu] = frmag(H,N);
z = poly(0,'z');
Hz=horner(H,1/z);
num = coeff(numer(Hz));
den = coeff(denom(Hz));
bStable = stable_EBSB(den(2:length(den)));
I = find(20*log10(abs(H_m)) < -3)
nu_c_est = nu(I(1));
endfunction
// -------------------------------------------------------------------------
function [bStable] = stable_EBSB(coeff_Den)
// Coeff_Den : contient les coefficients du dénominateur d'une fonction de transfert du type
// quotient de polynômes en z. Ne pas inclure le coefficient égal à 1
// les coefficients sont rangés suivant des puissances décroissantes.
// bStable : booléen
P = length(coeff_Den);
coeff2 = zeros(P,1);
for r = 1 : P
coeff2(r) = coeff_Den(P - r + 1);
end
mod_poles = abs(roots(poly([coeff2 ; 1],'z','coeff')));
I = find(mod_poles >= 1)
if(length(I) >= 1)
bStable = %F; // Pas stable au sens EBSB
else
bStable = %T; // Stable au sens EBSB
end
endfunction |
Partager