Высокочастотная фильтрация в MATLAB
Кто-нибудь знает, как использовать фильтры в MATLAB? Я не поклонник, поэтому меня не интересуют характеристики спада и т. Д. - у меня есть одномерный вектор сигнала x, сэмплированный на частоте 100 кГц, и я хочу выполнить фильтрацию верхних частот на нем (скажем, отклонить что-либо ниже 10 Гц) удалить базовый дрейф.
В справке описаны фильтры Баттерворта, Эллиптика и Чебычева, но нет простого объяснения того, как их реализовать.
3 ответа
Есть несколько фильтров, которые можно использовать, и фактический выбор фильтра будет зависеть от того, чего вы пытаетесь достичь. Поскольку вы упомянули фильтры Баттерворта, Чебышева и Эллиптика, я предполагаю, что вы ищете фильтры БИХ в целом.
Википедия - хорошее место, чтобы начать читать о различных фильтрах и о том, что они делают. Например, Баттерворт максимально плоский в полосе пропускания, а отклик в полосе останова скатывается. В Чебышеве у вас есть плавный отклик либо в полосе пропускания (тип 2), либо в полосе останова (тип 1), но и в более крупных нерегулярных рябях в другом, и, наконец, в эллиптических фильтрах имеются колебания в обеих полосах. Следующее изображение взято из Википедии.
Таким образом, во всех трех случаях вы должны обменять что-то на что-то другое. В Баттерворте у вас нет пульсаций, но частотная характеристика снижается медленнее. На приведенном выше рисунке, это берет от 0.4 до 0.55 чтобы получить половину власти. В Чебышеве вы получаете более крутой спад, но вы должны допускать неравномерную и более сильную рябь в одной из полос, а в Эллиптической вы получаете почти мгновенную отсечку, но у вас есть колебания в обеих полосах.
Выбор фильтра будет полностью зависеть от вашего приложения. Вы пытаетесь получить чистый сигнал практически без потерь? Тогда вам нужно что-то, что дает вам плавный отклик в полосе пропускания (Butterworth/Cheby2). Вы пытаетесь убить частоты в полосе задержания, и вы не будете возражать против незначительной потери в отклике в полосе пропускания? Тогда вам понадобится что-то гладкое в стоп-полосе (Cheby1). Нужны ли вам очень острые углы среза, т. Е. Что-нибудь, выходящее за пределы полосы пропускания, пагубно для вашего анализа? Если это так, вы должны использовать эллиптические фильтры.
О фильтрах БИХ следует помнить, что у них есть полюса. В отличие от КИХ-фильтров, где вы можете увеличить порядок фильтра с единственным последствием, являющимся задержкой фильтра, увеличение порядка БИХ-фильтров сделает фильтр нестабильным. Под нестабильным я имею в виду, что у вас будут полюсы, которые лежат за пределами круга юнитов Чтобы понять, почему это так, вы можете прочитать вики-статьи о фильтрах БИХ, особенно о стабильности.
Чтобы проиллюстрировать мою точку зрения, рассмотрим следующий полосовой фильтр.
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
РЕДАКТИРОВАТЬ
Это в ответ на комментарий @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 функция.
Дешевый способ сделать такую фильтрацию, которая не требует напряжения мозговых клеток в дизайне, нулях, полюсах, пульсации и тому подобном:
* 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-данные, что угодно.