Кто-нибудь знает, как использовать фильтры в MATLAB? Я не поклонник, поэтому меня не интересуют характеристики сползания и т.д. - У меня есть 1-мерный вектор сигнала x, отснятый на частоте 100 кГц, и я хочу выполнить фильтрацию высоких частот (например, отклонение чего-либо ниже 10 Гц), чтобы удалить базовый дрейф.
Существуют фильтры Butterworth, Elliptical и Chebychev, описанные в справке, но нет простого объяснения того, как реализовать.
Есть несколько фильтров, которые можно использовать, и фактический выбор фильтра будет зависеть от того, чего вы пытаетесь достичь. Поскольку вы упомянули фильтры Butterworth, Chebyschev и Elliptical, я предполагаю, что вы ищете фильтры IIR в целом.
Википедия - хорошее место, чтобы начать читать различные фильтры и что они делают. Например, Butterworth максимально плоский в полосе пропускания, и ответ откатывается в полосе остановки. В Chebyschev у вас есть гладкий ответ либо в полосе пропускания (тип 2), либо в полосе остановки (тип 1) и более крупной, нерегулярной ряби в другой и, наконец, в Эллиптические фильтры, рябь в обеих полосах. Следующее изображение взято из wikipedia.
Итак, во всех трех случаях вам нужно торговать чем-то другим. В Butterworth вы не получаете рябь, но частотная характеристика скатывается медленнее. На приведенном выше рисунке требуется от 0.4 до 0.55 до половины мощности. В Чебышеве вы становитесь более крутым, но вы должны допускать нерегулярные и большие ряби в одной из групп, а в Elliptical вы получаете почти мгновенный срез, но у вас есть рябь в обеих полосах.
Выбор фильтра будет полностью зависеть от вашего приложения. Вы пытаетесь получить чистый сигнал с небольшими потерями? Тогда вам нужно что-то, что дает вам плавный ответ в полосе пропускания (Butterworth/Cheby2). Вы пытаетесь убить частоты в полосе задержек, и вы не будете возражать против незначительной потери ответа в полосе пропускания? Тогда вам понадобится что-то гладкое в полосе остановки (Cheby1). Нужны ли вам чрезвычайно острые углы отсечения, т.е. Что-то немного за пределами полосы пропускания вредно для вашего анализа? Если это так, вы должны использовать эллиптические фильтры.
Что нужно помнить о фильтрах IIR, так это то, что у них есть полюса. В отличие от FIR-фильтров, где вы можете увеличить порядок фильтра с единственным разветвлением, являющимся задержкой фильтра, увеличение порядка фильтров IIR сделает фильтр неустойчивым. Нестабильно, я имею в виду, что у вас будут полюса, которые лежат вне единичного круга. Чтобы понять, почему это так, вы можете прочитать статьи wiki по фильтрам 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) для получения следующих
Чтобы получить более стабильный фильтр с вашими точными требованиями к дизайну, вам нужно будет использовать фильтры второго порядка, используя метод 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-объема, что угодно.