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

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

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

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

Википедия – хорошее место, чтобы начать читать по различным фильтрам и тем, что они делают. Например, Баттерворт максимально плоский в полосе пропускания, и ответ откатывается в полосе остановки. В Чебышеве у вас есть гладкий ответ в полосе пропускания (тип 2) или в полосе остановки (тип 1) и более крупной, нерегулярной ряби в другой и, наконец, в эллиптических фильтрах, рябь в обеих полосах. Следующее изображение взято из википедии .

введите описание изображения здесь

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

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

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

Чтобы проиллюстрировать мою точку зрения, рассмотрим следующий фильтр полосовых частот.

 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) чтобы получить следующее

введите описание изображения здесь

Чтобы получить более стабильный фильтр с вашими точными требованиями к дизайну, вам понадобится использовать фильтры второго порядка, используя метод zpk вместо ba , в 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 ansible на сайте обмена файлами filtfilthd

РЕДАКТИРОВАТЬ

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

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

введите описание изображения здесь

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

Создайте свой фильтр – например, используя [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 функцию filtfilt .

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

 * 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-изображения, данные 3D-объема, что угодно.

  • Запуск пользовательского приложения Android из браузера Android / Chrome
  • Что такое фильтры намерений в Android?
  • Почему изменение порядка стекирования зависит от зависания фильтра webkit?
  • Как использовать AutoCompleteTextView и заполнять его данными из веб-API?
  • Java8: сумма значений из определенного поля объектов в списке
  • Android - Внедрение поискового фильтра в RecyclerView
  • Как использовать фильтр объектов с поддержкой softlayer api?
  • ExtJs - Фильтрация сетки с полем поиска в заголовке столбца
  • Запуск пользовательского приложения для Android от браузера android
  • Пользовательский фильтр ArrayAdapter в ListView
  • Как эффективно фильтровать фрейм данных?
  • Давайте будем гением компьютера.