Progetto di filtri FIR a fase lineare con la tecnica MINMAX: listato del programma

In questo programma verrà sviluppata una tecnica di progetto di filtri FIR a fase lineare (in questo caso particolare passabasso e passabanda) utilizzando i metodi di realizzazione dei filtri ottimi che sfruttano l’approssimazione di Chebyshev e quindi il teorema dell’alternanza:

Teorema:

se P(jw) è una combinazione lineare di r funzioni coseno , cioè:

allora condizione necessaria e sufficiente perché P(jw) sia l’unica migliore approssimazione di una funzione continua D(jw) su A , sottoinsieme compatto di [0,p] , è che l’errore pesato :

E ( jw ) = W ( jw )[ D ( jw ) –P ( jw ) ]

presenti almeno r+1 frequenze di massimo o minimo in A , ovvero devono esistere in A (r+1) punti :

tali che: e tali che: e

OSS.: Il teorema permette di trovare la soluzione ottima secondo Chebyshev.

I filtri realizzati sono i seguenti:

─ basso.m : è un filtro passabasso di lunghezza pari a 13 punti e frequenze di transizione normalizzate rispetto alla frequenza di Nyquist pari a : f_p =0.333 e f_s =0.5 . La funzione peso invece moltiplica di un fattore 10 l’errore in banda passante;

─ banda.m : è un filtro passabanda di lunghezza pari a 21 punti e con bande attenuate comprese tra le frequenze normalizzate rispetto alla frequenza di Nyquist di : f₁ =0 , f₂ =0.2 e f₅ =0.8 , f₆ =1 e banda passante compresa tra f₃ =0.3 e f₄ =0.7 . La funzione peso invece moltiplica di un fattore 10 l’errore in banda passante;

─ bassopk.m : è un filtro passabasso di 61 punti e frequenze di transizione normalizzate pari a : f_p =0.2 e f_s =0.3. La funzione peso invece è unitaria in tutto il campo di frequenze.

La soluzione del problema secondo l’approssimazione di Chebyshev può essere ottenuta impostando l’algoritmo di Remez secondo Parks e McClellan come verrà proposto in questo programma, stando attenti che non sempre questo algoritmo converge.

Infine verrà visualizzata la risposta del filtro modificando la risoluzione in bit dei coefficienti del polinomio trigonometrico calcolati.

remez.m

% PROGETTO DI FILTRI FIR CON LA TECNICA MINMAX

% Verrà realizzato un filtro FIR a partire dalla

% sua funzione di trasferimento ideale introdotta

% dall'esterno considerando le SOLE frequenze positive

% il progetto tieni conto dell'algoritmo di Remez

% ed in particolare fa uso del solo caso 1 visto nella teoria ,in cui

% si considerano un numero di campioni pari ad NFILT=2r-1 (sempre dispari)

% ed il filtro è sempre a simmetria pari e deve essere composto da rettangoli

% in cui verranno specificati gli estremi di banda (filtro costante a pezzi)

% oppure deve essere un filtro lineare a pezzi

% r-1 sarà quindi l'ordine del polinomio trigonometrico mentre

% r+1 sono i massimi e minimi successivi dell'errore massimo di

% modulo uguale e segno alterno che individueremo (frequenze estremali)

% indico con EDGE il vettore contenente gli estremi di banda del filtro

% normalizzati alla frequenza di Nyquist (vettore con elementi compresi tra 0 e 1)

% WTX è il vettore contenente i pesi dell'errore normalizzati

% nelle bande specificate precedentemente in EDGE

% f è il vettore contenente le frequenze in cui conosciamo

% la funzione di trasferimento ideale di partenza

% EDGE,WTX,f devono essere introdotte dall'esterno tutte insieme con una function

% con la function 'peso' verrà individuato il valore della funzione peso,

% che si trova nel vettore WTX, in corrispondenza della banda del filtro

% in cui si trova l'elemento della griglia di frequenze preso in considerazione

% il polinomio trigonometrico nei punti della griglia si ottiene risolvendo

% un sistema di equazioni lineari che valuta i vari coefficienti del polinomio stesso

clear all

err=0;

while err==0

% i punti del filtro da introdurre devono corrispondere SOLO alle frequenze positive

F=input('inserisci il nome della function contenente il filtro ideale da realizzare: ','s');

if exist(F)==0

disp('LA FUNZIONE INTRODOTTA NON ESISTE');

else

err=1;

end

LGRID=input('introduci il valore della densità della griglia per calcolare E(w) : ');

grafici=input('Per disegnare i grafici ad ogni passo di interazione digitare "1" : ');

if grafici==1

DB=input('Per disegnare i grafici in dB digitare "1" : ');

else

DB=0;

end

[f,filtro,EDGE,WTX]=feval(F);% introduzione del filtro e delle sue caratteristiche

r=length(filtro);

strr=int2str(r);

disp(['Il numero di campioni del filtro da realizzare sono: ',strr]);

figure

hold on

plot(f,filtro,'r');

grafico=plot(f,filtro,'.b');

set(grafico,'markersize',14);

set(gcf,'color',[1 1 1]);

title('filtro ideale di partenza');

xlabel('f/(fc/2)');

zoom

NBANDS=length(WTX);

j=1; % costruzione della griglia di frequenze

for t=0:1/(LGRID*r):1

for i=1:2:NBANDS*2-1

if t>=EDGE(i) & t<=EDGE(i+1)

s(j)=t;

j=j+1;

end

L=length(s);

if s(1)>EDGE(1);

g(1)=EDGE(1);

g(2)=s(1);

k=3;

else

g(1)=s(1);

k=2;

end

indice1(1)=1;

kin=2;

for j=2:L-1 % aggiungiamo gli estremi di banda nella griglia

agg=0;

for i=1:2:NBANDS*2-1

if s(j)>EDGE(i) & s(j-1)<EDGE(i)

g(k)=EDGE(i);

indice1(kin)=k;

g(k+1)=s(j);

k=k+2;

kin=kin+1;

agg=1;

end

if s(j)<EDGE(i+1) & s(j+1)>EDGE(i+1)

g(k)=s(j);

g(k+1)=EDGE(i+1);

indice1(kin)=k+1;

k=k+2;

kin=kin+1;

agg=1;

end

if (s(j)==EDGE(i))|(s(j)==EDGE(i+1))

indice1(kin)=k;

kin=kin+1;

end

if agg==0

g(k)=s(j);

k=k+1;

end

g(k)=s(L);

indice1(kin)=k;

clear s

L=length(g);

indice2(1)=1; % costruzione della griglia nelle bande di transizione

kin=2;

if EDGE(1)>0

bt=(0:1/(LGRID*r):EDGE(1));

if bt(length(bt))<EDGE(1)

bt(length(bt)+1)=EDGE(1);

end

indice2(kin)=length(bt);

indice2(kin+1)=length(bt)+1;

kin=kin+2;

else

bt=[];

end

for j=1:NBANDS-1

bt=[bt,(EDGE(2*j):1/(LGRID*r):EDGE(2*j+1))];

if bt(length(bt))<EDGE(2*j+1)

bt(length(bt)+1)=EDGE(2*j+1);

end

indice2(kin)=length(bt);

indice2(kin+1)=length(bt)+1;

kin=kin+2;

end

if EDGE(2*NBANDS)<1

bt=[bt,(EDGE(2*NBANDS):1/(LGRID*r):1)];

if bt(length(bt))<1

bt(length(bt)+1)=1;

end

indice2(kin)=length(bt);

else

indice2=indice2(1:kin-2);

end

j=fix(L/r); % fissiamo a caso r+1 frequenze

wi(1)=g(1); % di partenza considerando la continua

for i=2:r+1

wi(i)=g((i-1)*j);

end

for j=2:r % inseriamo nelle frequenze estremali gli estremi di banda

for i=1:NBANDS

if wi(j)>EDGE(2*i-1) & wi(j-1)<EDGE(2*i-1)

wi(j)=EDGE(2*i-1);

elseif wi(j)<EDGE(2*i) & wi(j+1)>EDGE(2*i)

wi(j)=EDGE(2*i);

end

wi(r+1)=1;

passi=0;

cambiamenti=1; % comincia il ciclo di interazioni

while cambiamenti~=0 % che termina quando non ci sono più

cambiamenti=0; % cambiamenti nelle r+1 frequenze

passi=passi+1;

D=interlin(f,filtro,wi,r,r+1); % valori del filtro ideale nelle frequenze estremali

for j=1:r

COS(:,j)=(cos(wi*pi*(j-1)))';

end

for i=1:r+1

COS(i,r+1)=((-1)^(i-1))/peso(wi(i),EDGE,WTX,NBANDS);

end

% i valori che il filtro ideale assume nei punti della grilia si

% ottengono interpolando linearmente gli r punti di partenza :

% per far questo usiamo la function interlin.m

Y=COS\D; % risoluzione del sistema lineare

X=Y(1:r);

for i=1:L % valori del polinomio trigonometrico

p1(i)=X(1); % valutati nella griglia a partire

for j=2:r % dalle r+1 frequenze estremali

p1(i)=p1(i)+X(j)*cos(g(i)*pi*(j-1));

end

for i=1:r+1 % valori del polinomio trigonometrico

p2(i)=X(1); % nelle frequenze estremali

for j=2:r

p2(i)=p2(i)+X(j)*cos(wi(i)*pi*(j-1));

end

for i=1:length(bt) % valori del polinomio trigoniometrico

p3(i)=X(1); % nelle bande di transizione

for j=2:r

p3(i)=p3(i)+X(j)*cos(bt(i)*pi*(j-1));

end

% disegno della risposta in frequenza del filtro del filtro ad ogni passo se richiesto

strpassi=int2str(passi);

if grafici==1

figure

w1=p1;

w2=p2;

w3=p3;

if DB==1

for i=1:length(w1)

if abs(w1(i))<1e-5

w1(i)=1e-5;

end

for i=1:length(w2)

if abs(w2(i))<1e-5

w2(i)=1e-5;

end

for i=1:length(w3)

if abs(w3(i))<1e-5

w3(i)=1e-5;

end

w1=20*log10(abs(w1));

w2=20*log10(abs(w2));

w3=20*log10(abs(w3));

ylabel('dB');

end

hold on

for i=1:NBANDS

plot(g(indice1(2*i-1):indice1(2*i)),w1(indice1(2*i-1):indice1(2*i)),'k');

end

for i=1:length(indice2)/2

plot(bt(indice2(2*i-1):indice2(2*i)),w3(indice2(2*i-1):indice2(2*i)),'r');

end

grafico=plot(wi,w2,'.b');

set(grafico,'markersize',10);

set(gcf,'color',[1 1 1]);

xlabel('f/(fc/2)');

title(['risposta in frequenza del filtro al passo ',strpassi,' di interazione'])

axis([0 1 min([w1 w3]) max([w1 w3])])

zoom

hold off

end

for i=1:L % aggiornamento delle r+1 frequenze 'candidate'

% ad essere le r+1 frequenze estremali

w=peso(g(i),EDGE,WTX,NBANDS);

d=interlin(f,filtro,g(i),r,1);

E2=w*(d-p1(i));

w=peso(wi(1),EDGE,WTX,NBANDS);

d=interlin(f,filtro,wi(1),r,1);

E1=w*(d-p2(1));

w=peso(wi(r+1),EDGE,WTX,NBANDS);

d=interlin(f,filtro,wi(r+1),r,1);

E5=w*(d-p2(r+1));

if g(i)<wi(1)

if E2*E1>=0

if abs(E2)>abs(E1)

wi(1)=g(i);

p2(1)=p1(i);

cambiamenti=cambiamenti+1;

end

elseif abs(E2)>abs(E5)

for k=r:-1:1

wi(k+1)=wi(k);

p2(k+1)=p2(k);

end

wi(1)=g(i);

p2(1)=p1(i);

cambiamenti=cambiamenti+1;

end

j=1;

sc=0;

while j<=r & sc==0

if g(i)>=wi(j) & g(i)<wi(j+1)

w=peso(wi(j),EDGE,WTX,NBANDS);

d=interlin(f,filtro,wi(j),r,1);

E3=w*(d-p2(j));

w=peso(wi(j+1),EDGE,WTX,NBANDS);

d=interlin(f,filtro,wi(j+1),r,1);

E4=w*(d-p2(j+1));

if E2*E3>=0

if abs(E2)>abs(E3)

wi(j)=g(i);

p2(j)=p1(i);

sc=1;

cambiamenti=cambiamenti+1;

else

break

end

else

if abs(E2)>abs(E4)

wi(j+1)=g(i);

p2(j+1)=p1(i);

sc=1;

cambiamenti=cambiamenti+1;

else

break

end

else

j=j+1;

end

if g(i)>=wi(r+1)

if E2*E5>=0

if abs(E2)>abs(E5)

wi(r+1)=g(i);

p2(r+1)=p1(i);

cambiamenti=cambiamenti+1;

end

elseif abs(E2)>abs(E1)

for k=1:r

wi(k)=wi(k+1);

p2(k)=p2(k+1);

end

wi(r+1)=g(i);

p2(r+1)=p1(i);

cambiamenti=cambiamenti+1;

end

end % fine dell'aggiornamento delle r+1 frequenze estremali 'candidate'

strcambi=int2str(cambiamenti);

disp(['il numero di cambi effettuati al passo ',strpassi,' sono : ',strcambi]);

end % fine della ricerca delle frequenze estremali

% COSTRUZIONE DEI GRAFICI FINALI

if grafici==1 & DB==1

title('risposta in frequenza del filtro all''ultimo passo di iterazione in dB')

elseif grafici==1 & DB~=1

title('risposta in frequenza del filtro all''ultimo passo di iterazione')

end

w1=p1;

w2=p2;

w3=p3;

for i=1:length(w1)

if abs(w1(i))<1e-5

w1(i)=1e-5;

end

for i=1:length(w2)

if abs(w2(i))<1e-5

w2(i)=1e-5;

end

for i=1:length(w3)

if abs(w3(i))<1e-5

w3(i)=1e-5;

end

w1=20*log10(abs(w1));

w2=20*log10(abs(w2));

w3=20*log10(abs(w3));

if (grafici~=1) | (grafici==1 & DB==1)

figure

hold on

for i=1:NBANDS

plot(g(indice1(2*i-1):indice1(2*i)),p1(indice1(2*i-1):indice1(2*i)),'k');

end

for i=1:length(indice2)/2

plot(bt(indice2(2*i-1):indice2(2*i)),p3(indice2(2*i-1):indice2(2*i)),'r');

end

grafico=plot(wi,p2,'.b');

set(grafico,'markersize',10);

set(gcf,'color',[1 1 1]);

xlabel('f/(fc/2)');

title('risposta in frequenza del filtro all''ultimo passo di iterazione')

axis([0 1 min([p1 p3]) max([p1 p3])])

zoom

hold off

end

if (grafici~=1) | (grafici==1 & DB~=1)

figure

hold on

for i=1:NBANDS

plot(g(indice1(2*i-1):indice1(2*i)),w1(indice1(2*i-1):indice1(2*i)),'k');

end

for i=1:length(indice2)/2

plot(bt(indice2(2*i-1):indice2(2*i)),w3(indice2(2*i-1):indice2(2*i)),'r');

end

grafico=plot(wi,w2,'.b');

set(grafico,'markersize',10);

set(gcf,'color',[1 1 1]);

xlabel('f/(fc/2)');

ylabel('dB');

title('risposta in frequenza del filtro all''ultimo passo di iterazione in dB')

axis([0 1 min([w1 w3]) max([w1 w3])])

zoom

hold off

end

disp(['In totale i passi di iterazione sono: ',strpassi]);

strdelta=num2str(abs(Y(r+1)));

disp(['l''ampiezza massima delle oscillazioni normalizzata è : ',strdelta]);

% GRAFICI DELLA RISPOSTA IN FREQUENZA DEL FILTRO MODIFICANDO

% LA RISOLUZIONE DEI COEFFICIENTI DEL POLINOMIO TRIGONIOMETRICO

disp('se si desidera modificare la risoluzione dei coefficienti del polinomio');

risp=input('trigoniometrico digitare ''1'' : ');

if risp==1

figure

hold on

bit=4;

colori=['k','r','b'];

for ris=1:3

for i=1:r

xi(i)=floor((X(i)/max(X))*(2^bit-1))*(max(X)/(2^bit-1));

end

for i=1:L % valori del polinomio trigonometrico

p1(i)=xi(1); % valutati nella griglia a partire

for j=2:r % dalle r+1 frequenze estremali

p1(i)=p1(i)+xi(j)*cos(g(i)*pi*(j-1));

end

for i=1:length(bt) % valori del polinomio trigoniometrico

p3(i)=xi(1); % nelle bande di transizione

for j=2:r

p3(i)=p3(i)+xi(j)*cos(bt(i)*pi*(j-1));

end

for i=1:length(p1)

if abs(p1(i))<1e-5

p1(i)=1e-5;

end

for i=1:length(p3)

if abs(p3(i))<1e-5

p3(i)=1e-5;

end

p1=20*log10(abs(p1));

p3=20*log10(abs(p3));

if EDGE(1)>0

gr=bt(indice2(1):indice2(2));

w=p3(indice2(1):indice2(2));

j=3;

else

gr=[];

w=[];

j=1;

end

for i=0:NBANDS-2

gr=[gr,g(indice1(2*i+1):indice1(2*i+2)),bt(indice2(2*i+j):indice2(2*i+j+1))];

w=[w,p1(indice1(2*i+1):indice1(2*i+2)),p3(indice2(2*i+j):indice2(2*i+j+1))];

end

gr=[gr,g(indice1(2*NBANDS-1):indice1(2*NBANDS))];

w=[w,p1(indice1(2*NBANDS-1):indice1(2*NBANDS))];

if EDGE(2*NBANDS)<1

gr=[gr,bt(length(indice2)-1:length(indice2))];

w=[w,p3(indice2(length(indice2)-1):indice2(length(indice2)))];

end

plot(gr,w,colori(ris));

set(gcf,'color',[1 1 1]);

xlabel('f/(fc/2)');

ylabel('dB');

title(['risposta in frequenza del filtro con le seguenti risoluzioni'])

zoom

bit=bit*2;

end

legend('4 bit','8 bit','16 bit');

end

h(r)=X(1);

for i=2:r

h(r-i+1)=X(i)/2;

end

for j=r+1:2*r-1

h(j)=h(2*r-j);

end

disp('La risposta all''impulso del filtro è contenuta nel vettore "h" ');

h=h';