Filtro passa-alto in MATLAB

Qualcuno sa come usare i filtri in MATLAB? Non sono un aficionado, quindi non mi interessano le caratteristiche di roll-off, ecc. – Ho un vettore di segnale 1 x campionato a 100 kHz, e voglio eseguire un filtro passa-alto su di esso (ad esempio, rifiutando qualcosa sotto i 10Hz ) per rimuovere la deriva della linea di base.

Ci sono i filtri Butterworth, Ellittici e Chebychev descritti nell’aiuto, ma nessuna spiegazione semplice su come implementare.

Esistono diversi filtri che è ansible utilizzare e la scelta effettiva del filtro dipenderà da ciò che si sta tentando di ottenere. Dal momento che hai menzionato i filtri Butterworth, Chebyschev ed Ellittici, presumo che tu stia cercando filtri IIR in generale.

Wikipedia è un buon punto di partenza per leggere i diversi filtri e quello che fanno. Ad esempio, Butterworth è al massimo piatto nella banda passante e la risposta si interrompe nella banda di arresto. In Chebyschev , hai una risposta fluida sia nella banda passante (tipo 2) che nella banda di stop (tipo 1) e increspature più grandi e irregolari nell’altra e infine, nei filtri ellittici , ci sono increspature in entrambe le bande. L’immagine seguente è tratta da wikipedia .

inserisci la descrizione dell'immagine qui

Quindi in tutti e tre i casi, devi scambiare qualcosa per qualcos’altro. In Butterworth non ottieni increspature, ma la risposta in frequenza è più lenta. Nella figura sopra, ci vuole da 0.4 a circa 0.55 per arrivare a metà potenza. In Chebyschev, ottieni un rollover più ripido, ma devi tenere conto di increspature irregolari e più grandi in una delle bande, e in Ellittica, ottieni un taglio quasi istantaneo, ma hai delle increspature in entrambe le bande.

La scelta del filtro dipenderà interamente dalla tua applicazione. Stai cercando di ottenere un segnale pulito con perdite minime o nulle? Allora hai bisogno di qualcosa che ti dia una risposta fluida nella passabanda (Butterworth / Cheby2). Stai cercando di uccidere le frequenze nella banda di stop e non ti dispiacerà una piccola perdita nella risposta nella banda passante? Quindi avrai bisogno di qualcosa che sia fluido nella banda di stop (Cheby1). Avete bisogno di angoli di taglio estremamente nitidi, vale a dire, qualsiasi cosa al di là della banda passante è dannosa per la vostra analisi? Se è così, dovresti usare i filtri ellittici.

La cosa da ricordare sui filtri IIR è che hanno dei poli. A differenza dei filtri FIR in cui è ansible aumentare l’ordine del filtro con la sola ramificazione che è il ritardo del filtro, l’aumento dell’ordine dei filtri IIR renderà il filtro instabile. Con unstable, voglio dire che avrai pali che si trovano al di fuori del cerchio unitario. Per capire perché è così, puoi leggere gli articoli wiki sui filtri IIR , in particolare la parte sulla stabilità.

Per illustrare ulteriormente il mio punto, si consideri il seguente filtro passa-banda.

 fpass=[0.05 0.2];%# passband fstop=[0.045 0.205]; %# frequency where it rolls off to half power Rpass=1;%# max permissible ripples in stopband (dB) Astop=40;%# min 40dB attenuation n=cheb2ord(fpass,fstop,Rpass,Astop);%# calculate minimum filter order to achieve these design requirements [b,a]=cheby2(n,Astop,fstop); 

Ora se osservi il diagramma a polo zero usando zplane(b,a) , vedrai che ci sono diversi poli ( x ) che si trovano fuori dal cerchio dell’unità, il che rende instabile questo approccio.

inserisci la descrizione dell'immagine qui

e questo è evidente dal fatto che la risposta in frequenza è tutta in tilt. Usa freqz(b,a) per ottenere quanto segue

inserisci la descrizione dell'immagine qui

Per ottenere un filtro più stabile con i tuoi esatti requisiti di progettazione, dovrai utilizzare i filtri del secondo ordine utilizzando il metodo zpk invece di ba , in MATLAB. Ecco come per lo stesso filtro di cui sopra:

 [z,p,k]=cheby2(n,Astop,fstop); [s,g]=zp2sos(z,p,k);%# create second order sections Hd=dfilt.df2sos(s,g);%# create a dfilt object. 

Ora, se guardi alle caratteristiche di questo filtro, vedrai che tutti i poli si trovano all’interno del cerchio dell’unità (quindi sono stabili) e corrispondono ai requisiti di progettazione

inserisci la descrizione dell'immagine qui

inserisci la descrizione dell'immagine qui

L’approccio è simile per il butter e l’ ellip , con equivalenza di buttord e ellipord . La documentazione di MATLAB ha anche buoni esempi sulla progettazione di filtri. Puoi build su questi esempi e il mio per progettare un filtro in base a ciò che desideri.

Per utilizzare il filtro sui tuoi dati, puoi utilizzare filter(b,a,data) o filter(Hd,data) seconda del filtro che usi eventualmente. Se si desidera la distorsione di fase zero, utilizzare filtfilt . Tuttavia, questo non accetta oggetti dfilt . Quindi per il filtro a fase zero con Hd , utilizzare il file filtfilthd disponibile nel sito di scambio di file di Mathworks

MODIFICARE

Questo è in risposta al commento di @ DarenW. Il livellamento e il filtraggio sono due operazioni diverse, e sebbene siano simili per alcuni aspetti (la media mobile è un filtro passa-basso), non si può semplicemente sostituire l’uno con l’altro a meno che non si possa essere certi che non sarà preoccupazione nell’applicazione specifica.

Ad esempio, implementando il suggerimento di Daren su un segnale chirp lineare da 0-25 kHz, campionato a 100 kHz, questo spettro di frequenza dopo il livellamento con un filtro gaussiano

inserisci la descrizione dell'immagine qui

Certo, la deriva vicina a 10Hz è quasi nulla. Tuttavia, l’operazione ha completamente cambiato la natura delle componenti di frequenza nel segnale originale. Questa discrepanza si verifica perché hanno completamente ignorato il roll-off dell’operazione di livellamento (vedi linea rossa) e presumevano che sarebbe stato zero piatto. Se fosse vero, allora la sottrazione avrebbe funzionato. Ma ahimè, non è così, ed è per questo che esiste un intero campo sulla progettazione dei filtri.

Crea il tuo filtro – per esempio usando [B,A] = butter(N,Wn,'high') dove N è l’ordine del filtro – se non sei sicuro di cosa si tratta, basta impostarlo su 10. Wn è il cutoff frequenza normalizzata tra 0 e 1, con 1 corrispondente a metà della frequenza di campionamento del segnale. Se la frequenza di campionamento è fs e si desidera una frequenza di taglio di 10 Hz, è necessario impostare Wn = (10/(fs/2)) .

Puoi quindi applicare il filtro usando Y = filter(B,A,X) dove X è il tuo segnale. Puoi anche guardare nella funzione filtfilt .

Un modo cheapo per fare questo tipo di filtraggio che non coinvolge le cellule cerebrali tese su design, zeri e poli e increspature e tutto ciò, è:

 * Make a copy of the signal * Smooth it. For a 100KHz signal and wanting to eliminate about 10Hz on down, you'll need to smooth over about 10,000 points. Use a Gaussian smoother, or a box smoother maybe 1/2 that width twice, or whatever is handy. (A simple box smoother of total width 10,000 used once may produce unwanted edge effects) * Subtract the smoothed version from the original. Baseline drift will be gone. 

Se il segnale originale è spikey, potresti voler usare un filtro mediano corto prima del grande più liscio.

Ciò generalizza facilmente alle immagini 2D, ai dati del volume 3D, a prescindere.