[ Home page | Sommario | Ricerca | Invia | Rispondi | Avanti | Indietro | Su ]
Da: Giocosta
Remote User:
Date: Wednesday 20 December 2000
Time: 11.28.07
Questa è la mia versione dell'algoritmo di Remez, ovviamente rivolta a coloro che stanno studiando Elaborazione numerica dei segnali. Il programma è stato realizzato con matlab 5.3 % Progetto dei filtri FIR mediante l'algotitmo di Remez
% Definizione maschera N=input('numero (dispari) di campioni del filtro : '); fp=input('freq. di taglio passed-band (Hz): '); fs=input('freq. di taglio stopped-band (Hz): '); fc=input('freq. di campionamento (Hz): '); if fp>=fc/2 | fs>=fc/2 error('fp e fs devono essere maggiori della frequenza di Nyquist'); end
% Normalizzazione delle frequenze rispetto alla frequenza di campionamento fpnorm=fp/fc fsnorm=fs/fc
% Definizione del ripple in banda passante e in banda attenuata % N.B. Quello che si fissa è solo un rapporto relativo tra i due ripple e non % la loro entità in termini assoluti delta1=input('delta1 in banda passante= '); delta2=input('delta2 in banda attenuata= ');
nbit=input('quanti bit vuoi usare per il calcolo dei coefficienti? ');
L=(N-1)/2; r=L+1; % Suppongo di considerare una risoluzione in frequnza pari a 0.0001 freq=[0:0.0001:0.5]; risol=length(freq);
% Fisso le frequenze estremali iniziali metà in banda passante % e metà in banda attenuata. Fra di esse considero la continua, fp, fs e la freq. % di Nyquist, le altre saranno equispaziate. In tutto saranno r+1. clear Fes; fpnorm1=ceil(fpnorm*10000)/10000; k1=fpnorm1*10000; fsnorm1=ceil(fsnorm*10000)/10000; k2=fsnorm1*10000; Fes=[linspace(0,fpnorm1,ceil(L/2)+1) linspace(fsnorm1,0.5,floor(L/2)+1)];
% def. di funz. peso e funz. desiderata clear D; %Funzione desiderata clear W; %Funzione peso clear A; %Matrice dei coefficienti in A*bk=D clear H; W=[(delta2/delta1)*ones(1,length([0:0.0001:fpnorm1])-1) zeros(1,length(fpnorm1:0.0001:fsnorm1)-2) ones(1,length([fsnorm1:0.0001:0.5])+1)]; D=[ones(1,length([0:0.0001:fpnorm1])) zeros(1,length(fpnorm1:0.0001:0.5)-1)];
% Inizio dell'algoritmo vero e proprio
delta=0; iterate=0; condstop=0; while condstop==0 iterate=iterate+1; if iterate==26 error('il metodo sta richiedendo troppe iterate, quindi potrebbe non convergere'); end;
% def. matrice dei coefficienti a=length(find(Fes<=fpnorm1)); W1=[(delta2/delta1)*ones(1,a) ones(1,L+2-a)]; D1=[ones(1,a) zeros(1,L+2-a)]; for i=1:r+1 for j=1:r A(i,j)=cos(2*(j-1)*pi*Fes(i)); end; A(i,j+1)=(1/W1(i))*(-1)^(i); end;
% Calcolo della funzione di trasferimento del filtro % Risolvo prima il sistema bk=(inv(A)*D1'); bk=floor(bk/max(bk)*2^nbit)*max(bk)/2^nbit; H=zeros(1,risol); for n=1:r H=H+bk(n)*cos(2*pi*freq*(n-1)); end;
% Visualizzazione della funzione di trasferimento
figure plot(freq,H); hold on % Sovrapposizione di marcatori indicanti la posizione delle frequenze estremali kk=ceil(Fes*10000+1); plot(Fes,H(kk),'+'); hold off grid title('f.d.t.'); E=W.*(H-D);
%ricerca delle nuove freq. estremali (ricerca dei massimi e minimi locali) i=0; pos=[]; clear Fes1; for k=2:risol-1 %in questo modo calcolo i massimi e i minimi già ordinati if (H(k)>H(k-1) & H(k)>H(k+1)) | (H(k)<H(k-1) & H(k)<H(k+1)) i=i+1; %Bisogna tener presente sia la risoluzione e sia che %k non può partire da 0 come essendo un indice Fes1(i)=(k-1)*0.0001; end % Tengo conto della posizione delle frequenze estremali pos=[pos k]; end
% Fisso le nuove frequenze estremali
%Se il numero di max e min è r+1-4=r-3 vuol dire che devo aggiungere la continua, %fp,fc e fN, ma se ne mancano 3 allora devo scegliere se devo scegliere la %continua oppure fN. Se per esempio la continua e il primo max o min portano ad avere %un errore o positivo o negativo (per entrambi) allora vuol dire che molto probabilmente %la continua non sarà una freq estremale.Stessa cosa per fN e l'ultimo max o min trovato. %Se sia per la continua che per fN si hanno segni opposti (nel senso detto prima) %allora si sceglie tra le due quella che potrebbe dare origine ad uno scostamento %maggiore dalla maschera.
l=length(Fes1); %sarà anche la dimensione del vettore pos if l==r-3 %vale a dire r+1-4 Fes1=sort([Fes1 0 fpnorm1 fsnorm1 0.5]); % a questo punto devo scegliere tra la continua e Nyquist (sicuramente % saranno max e min locali, ora biosgna vedere se possono essere % frequnze estremali) elseif l==r-2 & sign(E(pos(1)))==sign(E(1)) & sign(E(pos(l)))~=sign(E(risol)) Fes1=sort([Fes1 fpnorm1 fsnorm1 0.5]); elseif l==r-2 & sign(E(pos(1)))~=sign(E(1)) & sign(E(pos(l)))==sign(E(risol)) Fes1=sort([Fes1 0 fpnorm1 fsnorm1 ]); elseif l==r-2 & abs(E(pos(1))-E(1))>abs(E(pos(l))-E(risol)) Fes1=sort([Fes1 0 fpnorm1 fsnorm1 ]); elseif l==r-2 & abs(E(pos(1))-E(1))<=abs(E(pos(l))-E(risol)) Fes1=sort([Fes1 fpnorm1 fsnorm1 0.5]);
else Fes1 error(' errore da qualche parte'); % Da togliere quando uscirà end
newdelta=bk(length(bk)); norma=abs(Fes-Fes1); norma2=max(norma);
if ((norma2<.0001) & abs(newdelta-delta)<.0001) condstop=1; end; Fes=Fes1; disp('Valore del delta alla iterata') iterata=iterate delta=newdelta end; %fine while esterno disp('frequenze estremali normalizzate a fc '); Fes disp('coefficienti della risposta all''impulso'); h=[bk(1:r)' bk(r-1:-1:1)'] iterate