Filtro passa-alto in MATLAB

Question

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 filtraggio passa-alto su di esso (ad esempio, rifiutando qualcosa sotto 10Hz) per rimuovere la deriva della linea di base.

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

Answer 1

Esistono diversi filtri che è possibile 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, si ha una risposta fluida sia nella banda passante (tipo 2) che nella banda di arresto (tipo 1) e increspature più grandi e irregolari nell'altra e, infine, nei filtri Elliptical, ci sono increspature in entrambe le bande. La seguente immagine is taken from wikipedia.

Quindi, in tutti e tre i casi, si deve scambiare qualcosa per qualcos'altro. In Butterworth non ottieni increspature, ma la risposta in frequenza è più lenta. Nella figura sopra, prende 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 dall'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, qualcosa che va oltre la banda passante è dannoso per la vostra analisi? Se è così, dovresti usare i filtri ellittici.

La cosa da ricordare sui filtri IIR è che hanno i poli. A differenza dei filtri FIR in cui è possibile 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 del wiki su IIR filters, in particolare la parte sulla stabilità.

Per illustrare ulteriormente il mio punto, considerare 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 si guarda il diagramma zero palo utilizzando zplane(b,a), vedrete che ci sono diversi poli (x) che si trovano al di fuori del cerchio unitario, il che rende questo approccio instabile.

e questo è evidente dal fatto che la risposta in frequenza è tutto in tilt. Utilizzare freqz(b,a) ottenere le seguenti

Per ottenere un filtro più stabile con i vostri requisiti esatti di progettazione, è necessario utilizzare i filtri del secondo ordine con il metodo z-p-k invece di b-a, in MATLAB. Ecco come per lo stesso filtro come 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 si guardano le caratteristiche di questo filtro, vedrete che tutti i poli si trovano all'interno del cerchio unitario (e quindi stabile) e corrisponde ai requisiti di progettazione

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

Per utilizzare il filtro sui dati, è possibile eseguire filter(b,a,data) o filter(Hd,data) a seconda del filtro che si utilizzerà. Se si desidera la distorsione di fase zero, utilizzare filtfilt. Tuttavia, questo non accetta gli oggetti dfilt. Quindi, per il filtro zero fase con Hd, utilizzare il file filtfilthd disponibili sul file sito di scambio Mathworks

EDIT

Questo è in risposta a @ 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, l'implementazione suggerimento di Daren su un segnale chirp lineare da 0-25kHz, campionato a 100kHz, questo lo spettro di frequenza dopo la rasatura con un filtro gaussiano

Certo, la deriva vicino 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.

Answer 2

Crea il tuo filtro - ad 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 è la frequenza di taglio 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)).

È quindi possibile applicare il filtro utilizzando Y = filter(B,A,X) dove X è il segnale. È anche possibile esaminare la funzione filtfilt.

Answer 3

Un modo cheapo per fare questo tipo di filtraggio che non coinvolge sforzare cellule del cervello sul design, zeri e poli e ondulazione e tutto ciò che è:

* 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, è potrebbe voler usare un filtro mediano corto prima del grande più liscio.

Questo generalizza facilmente a immagini 2D, dati di volume 3D, qualunque sia.