Высокочастотная фильтрация в MATLAB

Question

Кто-нибудь знает, как использовать фильтры в MATLAB? Я не поклонник, поэтому меня не интересуют характеристики сползания и т. Д. - У меня есть одномерный вектор сигнала x, отснятый на частоте 100 кГц, и я хочу выполнить фильтрацию высоких частот (например, отклонение чего-либо ниже 10 Гц), чтобы удалить базовый дрейф.

Существуют фильтры Butterworth, Elliptical и Chebychev, описанные в справке, но нет простого объяснения того, как реализовать.

Answer 1

Существует несколько фильтров, которые могут использоваться, и фактический выбор фильтра будет зависеть от того, чего вы пытаетесь достичь. Поскольку вы упомянули фильтры Butterworth, Chebyschev и Elliptical, я предполагаю, что вы ищете фильтры IIR в целом.

Википедия - это хорошее место для начала чтения различных фильтров и того, что они делают. Например, Butterworth максимально плоский в полосе пропускания, и ответ откатывается в полосе останова. В Chebyschev у вас есть гладкий отклик в полосе пропускания (тип 2) или в полосе остановки (тип 1) и большей, нерегулярной ряби в другой и, наконец, в фильтрах Elliptical, в обеих полосах пульсирует рябь. Следующее изображение is taken from wikipedia.

Так во всех трех случаях, вы должны что-то для чего-то еще торговать. В Butterworth вы не получаете рябь, но частотная характеристика скатывается медленнее. На приведенном выше рисунке требуется от 0.4 примерно до 0.55, чтобы добраться до половины мощности. В Чебышеве вы становитесь более крутым, но вы должны допускать нерегулярные и большие ряби в одной из групп, а в Elliptical вы получаете почти мгновенный срез, но у вас есть рябь в обеих полосах.

Выбор фильтра полностью зависит от вашего приложения. Вы пытаетесь получить чистый сигнал с небольшими потерями? Тогда вам нужно что-то, что дает вам плавный ответ в полосе пропускания (Butterworth/Cheby2). Вы пытаетесь убить частоты в полосе задержек, и вы не будете возражать против незначительной потери ответа в полосе пропускания? Тогда вам понадобится что-то гладкое в полосе остановки (Cheby1). Нужны ли вам чрезвычайно острые углы отсечения, т. Е. Что-то немного за пределами полосы пропускания вредно для вашего анализа? Если это так, вы должны использовать эллиптические фильтры.

Что нужно помнить о фильтрах IIR, так это то, что у них есть полюса. В отличие от FIR-фильтров, где вы можете увеличить порядок фильтра с единственным разветвлением, являющимся задержкой фильтра, увеличение порядка фильтров IIR сделает фильтр неустойчивым. Нестабильно, я имею в виду, что у вас будут полюса, которые лежат вне единичного круга. Чтобы понять, почему это так, вы можете прочитать статьи в вики на IIR filters, особенно в отношении стабильности.

Для дальнейшей иллюстрации рассмотрим следующий полосовой фильтр.

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); 

Теперь, если вы посмотрите на диаграмме нулевого полюса с помощью zplane(b,a), вы увидите, что есть несколько полюсов (x), лежащих вне единичного круга, что делает этот подход неустойчив.

и это видно из того, что частотная характеристика все наперекосяк. Используйте freqz(b,a), чтобы получить следующие

Чтобы получить более стабильный фильтр с вашими точными требованиями к конструкции, вы должны будете использовать фильтры второго порядка с использованием метода z-p-k вместо b-a в MATLAB. Вот как за тот же фильтр, как описано выше:

[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. 

Теперь, если вы посмотрите на характеристики этого фильтра, вы увидите, что все полюсы лежат внутри единичной окружности (отсюда стабильный) и соответствует требованиям к конструкции

подход аналогичен по butter и ellip, с эквивалентным buttord и ellipord. В документации MATLAB также есть хорошие примеры для проектирования фильтров. Вы можете использовать эти примеры и мои, чтобы создать фильтр в соответствии с тем, что вы хотите.

Чтобы использовать фильтр для ваших данных, вы можете сделать filter(b,a,data) или filter(Hd,data) в зависимости от того, какой фильтр вы в конечном итоге используете. Если вы хотите искажение нулевой фазы, используйте filtfilt. Однако это не принимает dfilt объектов. Таким образом, чтобы нулевой фазовый фильтр с Hd, используйте filtfilthd файл доступный на файл обмена сайта Mathworks

EDIT

Это в ответ на @ комментарий DarenW в. Сглаживание и фильтрация - это две разные операции, и хотя они схожи в некоторых отношениях (скользящее среднее - это фильтр нижних частот), вы не можете просто заменить один на другой, если только вы не можете быть уверены, что он не будет в конкретном приложении.

Например, реализация предложения Дарен на линейный сигнал с ЛЧМ от 0-25kHz, дискретизированного с частотой 100 кГц, этот частотный спектр после сглаживания с гауссовым фильтром

Конечно, дрейф близок к 10Гц почти равна нулю. Однако операция полностью изменила характер частотных составляющих исходного сигнала. Это расхождение происходит потому, что они полностью игнорировали спуск операции сглаживания (см. Красную линию) и предполагали, что он будет плоским нолем. Если бы это было так, то вычитание сработало бы. Но, увы, это не так, поэтому существует целая область разработки фильтров.

Answer 2

Создайте свой фильтр - например, используя [B,A] = butter(N,Wn,'high'), где N - порядок фильтра - если вы не знаете, что это такое, просто установите его на 10. Wn - частота среза, нормированная между 0 и 1, с 1, соответствующей половину частоты дискретизации сигнала. Если ваша частота дискретизации равна fs, и вам нужна частота среза 10 Гц, вам необходимо установить Wn = (10/(fs/2)).

Затем вы можете применить фильтр, используя Y = filter(B,A,X), где X - ваш сигнал. Вы также можете посмотреть функцию filtfilt.

Answer 3

дешевки способ сделать этот вид фильтрации, которая не предполагает напрягает клетки мозга на дизайн, нулей и полюсов и пульсация, и все, что есть:

* 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. 

Если исходный сигнал Spikey, вы возможно, захотите использовать короткий срединный фильтр перед большей плавностью.

Это легко обобщает 2D-изображения, данные трехмерного объема, что угодно.