Filtro digitale prima prova

[DirtyDeeds]

e una frequenza di campionamento

Volevi dire ?

Se è così non va bene, devi scegliere .

[simo85]

Si ho sbagliato a scrivere.

Grazie per la precisazione!

[DirtyDeeds]

Tra l'altro c'è qualcosa che non mi quadra nelle equazioni per il calcolo dei coefficienti che hai scritto: perché quelle siano valide direi che dovresti avere:

[dimaios]

simo85, a parte il problema delle frequenze normalizzate dovrebbe insospettirti il fatto che la risposta impulsiva che hai ottenuto non è simmetrica.
Ho provato a progettare un FIR di ordine 32 con frequenza normalizzata pari a 0.5 ed il risultato è il seguente.

CoefficientsFIR.jpg (20.09 KiB) Osservato 2114 volte

I coefficienti devono essere simmetrici e tutti nella parte altrimenti il filtro sarebbe anticausale ovvero richiederebbe i campioni futuri.
Questa cosa non è impossibile nel senso che se hai già i tuoi campioni memorizzati da qualche parte e vuoi filtrarli off-line è possibile avere a disposizione i campioni successivi a quello attuale.
Nel tuo caso invece vuoi elaborare il segnale d'ingresso in realtime per cui la risposta impulsiva deve essere causale.

Ad ogni modo una frequenza di campionamento pari a 10 kHz implica che il range di frequenze "filtrabili" digitalmente è tra 0 e 5 kHz.
Se scegli il taglio proprio a 5 kHz vuol dire che ottieni di fatto un passa-tutto.

Siccome il campionamento di un segnale in frequenza si traduce nella ripetizione periodica del suo spettro, il filtro continuo che vuoi discretizzare ( oppure direttamente sintetizzare nel dominio discreto ) potrà essere specificato da e .
In pratica la frequenza di taglio del tuo passa-basso dovrà essere compresa tra gli estemi precedentemente citati.
Se ottieni di fatto un filtro passa-tutto.

Sorgente FidoCadJ | Ingrandisci | Cos'è Fidocad

Per fare l'esercizio imposta ed .
In questo modo potrai eseguire un test di validazione dando in ingresso al filtro la somma di due sinusoidi, una in banda passante e l'altra al di fuori di essa.



ed .

Se tutto funziona correttamente all'uscita dovresti ottenere un segnale che approssima molto bene il segnale .

[dimaios]

simo85 ti invito ad impiegare un tool per la sintesi del filtro così prendi mano con i parametri di design e solo successivamente andrai ad indagare le tecniche di progettazione.

Per esempio se usi questo hai in uscita già i coefficienti in double oppure in int ... anzi ... fornisce addirittura il codice già scritto che puoi facilmente studiare.

Il vantaggio è quello di concentrarsi sull'implementazione ed i relativi problemi.
Poi puoi indagare a fondo sulle tecniche di sintesi.
Ce ne sono moltissime, ognuna delle quali ha pro e contro che devono essere valutati in base all'applicazione.

Anche se non banale la miglior tecnica di progettazione per i filtri FIR è sicuramente quella di Parks McClellan che utilizza l'algoritmo di Remez.

[DirtyDeeds]

In pratica la frequenza di taglio del tuo passa-basso dovrà essere compresa tra gli estemi precedentemente citati.
Se ottieni di fatto un filtro passa-tutto.

Infatti. Con le notazioni di simo85, per , si deve allora avere



cosa che non viene se si definisce



Con la pulsazione



invece, per , si ha e



e



per .

Per evitare questo tipo di pasticci, e per evitare di avere che gira nei coefficienti, nell'analisi dei segnali preferisco sempre lavorare con le frequenze e non con le pulsazioni.

[simo85]

Urca, quindi il libro su cui ho fatto riferimento mi ha dato una cantonata?

Menomale che l'ho trovato "in biblioteca" (per così dire)..

Ora sono appena tornato a casa, devo fare un paio di cose in casa, successivamente vedo di aggiustare il tutto secondo i suggerimenti dati..

Vada per la tool offline ovviamente. Dò anche una occhiata al sito linkato, e vedo se riesco a fare il tutto anche con Octave.

Però una cosa la voglio dire. Mi sono fatto 200 pagine di libro studiate di volta in volta in metropolitana la mattina dopo il caffè, e più va avanti questo thread che è appena cominciato, più mi torna in mente l'importanza della pratica assieme alla teoria. Ma questo in tutto.

Siete fantastici, per me è praticamente nuovo l'argomento e passo dopo passo ci metto il mio impegno. Dovevo dirlo..

Vedo di mettermi al lavoro.

[DirtyDeeds]

Urca, quindi il libro su cui ho fatto riferimento mi ha dato una cantonata?

Gli errori di stumpa sono sempre in agguato Caveat emptor: prima di usare delle equazioni per un lavoro, da qualunque fonte o autore esse arrivino, se proprio non le si vuole riderivare, fare sempre un controllo su casi semplici per verificare che diano risultati consistenti, di errori ne facciamo tutti

Tieni presente che un segnale discreto non è nient'altro che una successione di numeri, il tempo è sostituito da un indice intero: nel dominio discreto, la frequenza è una grandezza adimensionata definita nell'intervallo . Come nel dominio continuo, la pulsazione è data da .

L'analisi di Fourier, continua o discreta, si può fare sia nel dominio delle frequenze che in quello delle pulsazioni, in quest'ultimo caso portandosi dietro (a fatica) un po' di coefficienti . I matematici spesso preferiscono il dominio delle pulsazioni, ma tra gli ingegneri direi che prevale il dominio delle frequenze.

[simo85]

Eccomi.
Dunque, se



e la lunghezza dell'array h (la chiamerò N_h) deve essere la metà di quella dell'array x e y (ma forse non necessariamente - non mi è ancora del tutto chiaro - ) e

nel dominio discreto, la frequenza f è una grandezza adimensionata definita nell'intervallo [-1/2,1/2].

Rappresento la risposta all' impulso h per N = 32 punti ed N_h = N/2.

Sorgente FidoCadJ | Ingrandisci | Cos'è Fidocad

A seguire ho scritto il programma in C per stampare i coefficienti secondo la risposta d'impulso.
Ho pure corretto la condizione di input del numero all'esecuzione del programma (prima dava Segmentation Fault per pure questioni di endianess ed allocazioni di memoria).
Ho anche fatto in modo che all'esecuzione sia possibile scegliere individualmente la lunghezza di x[] e la lunghezza di h[]

Per fare l'esercizio imposta F_c=10kHz ed F_t=1kHz.

OK quindi e .

Codice: Seleziona tutto

    #include <math.h>
    #include <stdio.h>
    #include <stdlib.h>

    /* Pi costant*/
    #define PI 3.1415926535897932384626433832795029L

    void zeros(double x[], double h[], double y[], int N, int N_h)
    {
        unsigned int i = 0, tmpN = 0;
       
        if(N > N_h)
        {
            while(i < N)
            {   
                if(i < N_h)
                    h[i] = 0;

                x[i] = y[i] = 0;
           
                i++;
            }
        }     
        else
        {
            while(i < N_h)
            {   
                if(i < N)
                    x[i] = y[i] = 0;
                       
                h[i] = 0;

                i++;
            }
        }
    }

    double set_LP_coefficient(int n_i, double w_c)
    {
        if(n_i == 0)
            return 2 * w_c;
        else
            return sin(2 * PI * w_c * n_i)/(PI * n_i);
    }

    int main(int argc, char *argv[])
    {   
        if(argc < 2)
        {
            printf("usage: %s N(int)  N_h (int)\n\non is the \
    lenght of x array\nN_h is the lenght of array h\n\n", argv[0]);
            return 1  ;
        }
        else
        {
            /*
                i is the index for arrays
                n_i is the n point
                N is the total point number
                N_h is the impulse response lenght
             */
            int i = 0, n_i = 0;   
            int N = atoi(argv[1]), N_h = atoi(argv[2]);

            /*
                x[N] is the input signal array
                h[N_h] is the impulse response coefficient array
                y[N] is the output signal array

            */
            double x[N], h[N_h], y[N];

            /*
                fs Sampling frecuency set to 10kHz
                fc Cutoff frequency set to 1kHz
                w is 2PI(fc / fs)
            */
            double fc = 1000, fs = 10000, w_c = ((2 * PI * fc) / fs);

            zeros(x, h, y, N, N_h);

            /* set n_i to negative and N_h / 2) + 1 */
            n_i = 0 - (N_h / 2);

            /* set index to zero */
            i = 0;

            while(i < N_h)
            {
                h[i] = set_LP_coefficient(n_i, w_c);
                printf("h[%.2d]\t%.10f\n", i, h[i]);
               
                n_i++;
                i++;
            }
           
            return 0;
        }
    }



A seguire l'esecuzione del programma con N = 32 ed N_h = 16

Codice: Seleziona tutto

h[00]   0.0431684235
h[01]   -0.0350947863
h[02]   -0.0273938263
h[03]   0.0795082871
h[04]   -0.0375186924
h[05]   -0.1148614620
h[06]   0.2927940834
h[07]   0.2000000000
h[08]   0.2927940834
h[09]   -0.1148614620
h[10]   -0.0375186924
h[11]   0.0795082871
h[12]   -0.0273938263
h[13]   -0.0350947863
h[14]   0.0431684235
h[15]   -0.0033146847


Ora direi che è simmetrica.

Spero non ci siano errori.

Domani vedo di suddividere il programma in funzioni, e se è tutto OK vedo di procedere alla stesura del codice FIR, ora però vado a nanna che è tardi, ho i neuroni in tilt ed alle sette mi alzo.

EDIT: in rete ho trovato questo che sembra essere altrettanto interessante.

[dimaios]

e la lunghezza dell'array h (la chiamerò N_h) deve essere la metà di quella dell'array x e y (ma forse non necessariamente - non mi è ancora del tutto chiaro - ) e

simo85, la lunghezza della sequenza non dipende dal numero di campioni in ingresso.

Se per esempio scegli di fare un semplice filtro di media con 2 campioni passati avrai :



Questo e' un filtro FIR ( infatti la risposta all'impulso e' finita ).
In particolare la risposta impulsiva vale 0.5 negli istanti temporali 0 e 1 mentre risulta nulla altrove.
Come vedi con solo 2 coefficienti hai un filtro FIR che esegue la media di due campioni.

In generale il numero di coefficienti del filtro deriva dalle specifiche in frequenza del medesimo.
Se vuoi ottenere un filtraggio piu' spinto il numero di coefficienti aumentera'.

Prima di lanciare la simulazione in C ti conviene vedere lo spettro del filtro per avere un'idea dell'attenuazione e fare una previsione su cosa aspettarti in uscita ( magari utilizza octave ).