В этом примере показано, как спроектировать, анализируйте и примените цифровой фильтр к своим данным. Это поможет вам ответить на вопросы, такие как: как я компенсирую задержку, введенную фильтром?, Как я стараюсь не искажать свой сигнал?, Как я удаляю нежелательное содержимое из своего сигнала?, Как я дифференцирую свой сигнал?, и Как я интегрирую свой сигнал?
Фильтры могут использоваться, чтобы сформировать спектр сигнала желаемым способом или выполнить математические операции, такие как дифференцирование и интегрирование. В дальнейшем вы изучите некоторые практические концепции, которые упростят использование фильтров, когда вам будут нужны они.
Этот пример фокусируется на приложениях цифровых фильтров, а не на их проекте. Если вы хотите узнать больше, как спроектировать цифровые фильтры, смотрите Практическое Введение в пример Создания цифровых фильтров.
Цифровые фильтры вводят задержку вашего сигнала. В зависимости от характеристик фильтра задержка может быть постоянной по всем частотам, или она может меняться в зависимости от частоты. Тип задержки определяет меры, которые необходимо принять, чтобы компенсировать его. Функция grpdelay позволяет вам смотреть на задержку фильтра как на функцию частоты. Рассмотрение выхода этой функции позволяет вам идентифицировать, является ли задержка фильтра постоянной или если это меняется в зависимости от частоты (т.е. если это зависимо частотой).
Отфильтруйте задержку, которая является постоянной по всем частотам, может быть легко компенсирован путем сдвига сигнала вовремя. КИХ-фильтры обычно имеют постоянную задержку. С другой стороны, задержитесь, который меняется в зависимости от частоты, вызывает искажение фазы и может значительно изменить форму волны сигнала. Компенсация зависимой частотой задержки не так тривиальна что касается постоянного случая задержки. БИХ-фильтры вводят зависимую частотой задержку.
Компенсация постоянной задержки фильтра
Как упомянуто прежде, можно измерить группу задержки фильтра, чтобы проверить, что это - постоянная функция частоты. Можно использовать функцию grpdelay, чтобы измерить задержку фильтра, D, и компенсировать эту задержку путем добавления D нулей к входному сигналу и сдвига выходного сигнала вовремя выборками D.
Полагайте, что шумная электрокардиограмма сигнализирует, что вы хотите отфильтровать, чтобы удалить высокочастотный шум выше 75 Гц. Вы хотите применить КИХ фильтр lowpass и компенсировать задержку фильтра так, чтобы шумные и отфильтрованные сигналы были выровнены правильно и могли быть построены друг на друге для сравнения.
Fs = 500; % sample rate in Hz N = 500; % number of signal samples rng default; x = ecg(N)'+0.25*randn(N,1); % noisy waveform t = (0:N-1)/Fs; % time vector % Design a 70th order lowpass FIR filter with cutoff frequency of 75 Hz. Fnorm = 75/(Fs/2); % Normalized frequency df = designfilt('lowpassfir','FilterOrder',70,'CutoffFrequency',Fnorm);
Постройте групповую задержку фильтра, чтобы проверить, что это постоянно через все частоты, указывающие, что фильтр является линейной фазой. Используйте групповую задержку, чтобы измерить задержку фильтра.
grpdelay(df,2048,Fs) % plot group delay D = mean(grpdelay(df)) % filter delay in samples
D =
35
Перед фильтрацией добавьте нули D в конце вектора входных данных, x. Это гарантирует, что все полезные выборки сбрасываются из фильтра, и что входной сигнал и компенсированный задержке выходной сигнал имеют ту же длину. Отфильтруйте данные и компенсируйте задержку путем сдвига выходного сигнала выборками D. Этот последний шаг эффективно удаляет переходный процесс фильтра.
y = filter(df,[x; zeros(D,1)]); % Append D zeros to the input data y = y(D+1:end); % Shift data to compensate for delay figure plot(t,x,t,y,'r','linewidth',1.5); title('Filtered Waveforms'); xlabel('Time (s)') legend('Original Noisy Signal','Filtered Signal'); grid on axis tight
Компенсация зависимой частотой задержки
Зависимая частотой задержка вызывает искажение фазы в сигнале. Компенсация этого типа задержки не так тривиальна что касается постоянного случая задержки. Если ваше приложение позволяет оффлайн обрабатывать, можно удалить зависимую частотой задержку путем реализации фильтрации нулевой фазы с помощью функции filtfilt. filtfilt выполняет фильтрацию нулевой фазы путем обработки входных данных и в прямых и в обратных направлениях. Основной эффект состоит в том, что вы получаете искажение нулевой фазы, т.е. вы фильтруете данные с эквивалентным фильтром, который имеет постоянную задержку 0 выборок. Другие эффекты состоят в том, что вы получаете передаточную функцию фильтра, которая равняется величине в квадрате исходной передаточной функции фильтра и порядку фильтра, который удваивает порядок исходного фильтра.
Считайте сигнал ECG заданным в предыдущем разделе. Отфильтруйте этот сигнал с и незамедлительно компенсацию.
% Design a 7th order lowpass IIR elliptic filter with cutoff frequency % of 75 Hz. Fnorm = 75/(Fs/2); % Normalized frequency df = designfilt('lowpassiir',... 'PassbandFrequency',Fnorm,... 'FilterOrder',7,... 'PassbandRipple',1,... 'StopbandAttenuation',60);
Постройте групповую задержку фильтра и заметьте, что это меняется в зависимости от частоты, указывающей, что задержка фильтра зависима частотой.
grpdelay(df,2048,'half',Fs)
Отфильтруйте данные и посмотрите на эффекты каждой реализации фильтра на сигнале времени.
y1 = filter(df,x); % non-linear phase filter - no delay compensation y2 = filtfilt(df,x); % zero-phase implementation - delay compensation figure plot(t,x); hold on plot(t,y1,'r','linewidth',1.5); plot(t,y2,'linewidth',1.5); title('Filtered Waveforms'); xlabel('Time (s)') legend('Original Signal','Non-linear phase IIR output',... 'Zero-phase IIR output'); ax = axis; axis([0.25 0.55 ax(3:4)]) grid on
Заметьте, как нулевая фаза, фильтрующая эффективно, удаляет задержку фильтра.
Фильтрация нулевой фазы является большим инструментом, если ваше приложение допускает непричинные прямые/обратные операции фильтрации, и для изменения ответа фильтра на квадрат исходного ответа.
Фильтры, которые вводят постоянную задержку, являются линейными фильтрами фазы. Фильтры, которые вводят зависимую частотой задержку, являются нелинейными фильтрами фазы.
Фильтры обычно используются, чтобы удалить нежелательное спектральное содержимое из сигнала. Можно принять решение из множества фильтров сделать это. Вы выбираете фильтр lowpass, когда это необходимо, чтобы удалить высокочастотное содержимое или фильтр highpass, когда это необходимо, чтобы удалить низкочастотное содержимое. Можно также выбрать полосовой фильтр, чтобы удалить низкую частоту и высокочастотное содержимое при оставлении промежуточной полосы частот неповрежденной. Вы выбираете заграждающий фильтр, когда это необходимо, чтобы удалить частоты по данной полосе.
Рассмотрите звуковой сигнал, который имеет гул линии электропередачи и белый шум. Гул линии электропередачи вызывается тоном на 60 Гц. Белый шум является сигналом, который существует через всю полосу частот аудиосигнала.
Загрузите звуковой сигнал.
Fs = 44100; % Sample rate y = audioread('noisymusic.wav');
Постройте спектр мощности сигнала. Красный треугольный маркер показывает сильный тон на 60 Гц, вмешивающийся в звуковой сигнал.
[P,F] = pwelch(y,ones(8192,1),8192/2,8192,Fs,'power'); helperFilterIntroductionPlot1(F,P,[60 60],[-9.365 -9.365],... {'Original signal power spectrum', '60 Hz Tone'})
Можно сначала удалить как можно больше белого шумового спектрального содержимого использование фильтра lowpass. Полоса пропускания фильтра должна быть установлена в значение, которое предлагает хороший компромисс между шумоподавлением и аудио ухудшением из-за потери высокочастотного содержимого. Применение фильтра lowpass прежде, чем удалить гул на 60 Гц очень удобно, поскольку вы сможете проредить ограниченный полосой сигнал. Более низкий сигнал уровня позволит вам проектировать более резкий и более узкий заграждающий фильтр на 60 Гц с меньшим порядком фильтра.
Спроектируйте фильтр lowpass с частотой полосы пропускания 1 кГц и частотой полосы задерживания 1,4 кГц. Выберите проект минимального порядка.
Fp = 1e3; % Passband frequency in Hz Fst = 1.4e3; % Stopband frequency in Hz Ap = 1; % Passband ripple in dB Ast = 95; % Stopband attenuation in dB % Design the filter df = designfilt('lowpassfir','PassbandFrequency',Fp,... 'StopbandFrequency',Fst,'PassbandRipple',Ap,... 'StopbandAttenuation',Ast,'SampleRate',Fs); % Analyze the filter response hfvt = fvtool(df,'Fs',Fs,'FrequencyScale','log',... 'FrequencyRange','Specify freq. vector','FrequencyVector',F);
% Filter the data and compensate for delay D = mean(grpdelay(df)); % filter delay ylp = filter(df,[y; zeros(D,1)]); ylp = ylp(D+1:end); close(hfvt)
Посмотрите на спектр lowpass отфильтрованный сигнал. Отметьте, как содержимое частоты выше 1 400 Гц было удалено.
[Plp,Flp] = pwelch(ylp,ones(8192,1),8192/2,8192,Fs,'power'); helperFilterIntroductionPlot1(F,P,Flp,Plp,... {'Original signal','Lowpass filtered signal'})
Из графика спектра мощности выше, вы видите, что максимальное ненезначительное содержимое частоты lowpass отфильтрованный сигнал на уровне 1 400 Гц. Теоремой отсчетов демонстрационной частотой 2*1400 = 2 800 Гц были бы достаточны, чтобы представлять сигнал правильно, вы однако, используют частоту дискретизации 44 100 Гц, которая является отходами, поскольку необходимо будет обработать больше выборок, чем необходимые. Можно проредить сигнал уменьшать частоту дискретизации и уменьшать вычислительную загрузку путем сокращения количества выборок, которые необходимо обработать. Более низкая частота дискретизации также позволит вам проектировать более резкий и более узкий заграждающий фильтр, должен был удалить шум на 60 Гц, с меньшим порядком фильтра.
Downsample lowpass отфильтрованный сигнал фактором 10, чтобы получить частоту дискретизации Фс/10 = 4,41 кГц. Постройте спектр сигнала до и после субдискретизации.
Fs = Fs/10; yds = downsample(ylp,10); [Pds,Fds] = pwelch(yds,ones(8192,1),8192/2,8192,Fs,'power'); helperFilterIntroductionPlot1(F,P,Fds,Pds,... {'Signal sampled at 44100 Hz', 'Downsampled signal, Fs = 4410 Hz'})
Теперь удалите тон на 60 Гц с помощью БИХ-заграждающего фильтра. Позвольте полосе задерживания иметь ширину 4 Гц, сосредоточенных на уровне 60 Гц. Мы выбираем фильтр IIR, чтобы достигнуть резкой метки частоты, маленькой неравномерности в полосе пропускания и относительно младшего разряда. Обработайте данные с помощью filtfilt, чтобы избежать искажения фазы.
% Design the filter df = designfilt('bandstopiir','PassbandFrequency1',55,... 'StopbandFrequency1',58,'StopbandFrequency2',62,... 'PassbandFrequency2',65,'PassbandRipple1',1,... 'StopbandAttenuation',60,'PassbandRipple2',1,... 'SampleRate',Fs,'DesignMethod','ellip'); % Analyze the magnitude response hfvt = fvtool(df,'Fs',Fs,'FrequencyScale','log',... 'FrequencyRange','Specify freq. vector','FrequencyVector',Fds(Fds>F(2)));
Выполните фильтрацию нулевой фазы, чтобы избежать искажения.
ybs = filtfilt(df,yds);
Наконец, сверхдискретизируйте сигнал возвратить его исходной частоте дискретизации звука 44,1 кГц, которая совместима с аудио звуковыми картами.
yf = interp(ybs,10); Fs = Fs*10;
Бросьте итоговый взгляд на спектр исходных и обработанных сигналов. Заметьте, как высокочастотный уровень шума и тон на 60 Гц были ослаблены фильтрами.
[Pfinal,Ffinal] = pwelch(yf,ones(8192,1),8192/2,8192,Fs,'power'); close(hfvt) helperFilterIntroductionPlot1(F,P,Ffinal,Pfinal,... {'Original signal','Final filtered signal'})
Слушайте сигнал до и после обработки. Как упомянуто выше, конечный результат - то, что вы эффективно ослабили гул на 60 Гц и высокочастотный шум на звуковом файле.
% Play the original signal hplayer = audioplayer(y, Fs); play(hplayer); % Play the noise-reduced signal hplayer = audioplayer(yf, Fs); play(hplayer);
Функция diff MATLAB дифференцирует сигнал с недостатком, что можно потенциально увеличить уровень шума при выходе. Более оптимальный вариант состоит в том, чтобы использовать фильтр дифференциатора, который действует как дифференциатор в полосе интереса, и как аттенюатор на всех других частотах, эффективно удаляя высокочастотный шум.
Как пример, анализируйте скорость смещения пола создания во время землетрясения. Смещение или измерения дрейфа были зарегистрированы на первом этаже трех тестовых структур истории при условиях землетрясения и сохраненные в quakedrift.mat файле. Длина вектора данных 10e3, частота дискретизации составляет 1 кГц, и модули измерений являются cm.
Дифференцируйте данные о смещении, чтобы получить оценки скорости и ускорение пола создания во время землетрясения. Сравните результаты с помощью diff и КИХ-фильтра дифференциатора.
load quakedrift.mat Fs = 1000; % sample rate dt = 1/Fs; % time differential t = (0:length(drift)-1)*dt; % time vector
Спроектируйте 50-й фильтр дифференциатора порядка с частотой полосы пропускания 100 Гц, которая является пропускной способностью, по которой найдена большая часть энергии сигнала. Установите частоту полосы задерживания фильтра к 120 Гц.
df = designfilt('differentiatorfir','FilterOrder',50,... 'PassbandFrequency',100,'StopbandFrequency',120,... 'SampleRate',Fs);
Различная функция может рассматриваться как КИХ-фильтр первого порядка с ответом
. Используйте FVTool, чтобы сравнить ответ величины 50-го КИХ-фильтра дифференциатора порядка и ответ различной функции. Безусловно, оба ответа эквивалентны в области полосы пропускания (от 0 до 100 Гц). Однако в области полосы задерживания, 50-й фильтр порядка ослабляет компоненты, в то время как различный ответ усиливает компоненты. Это эффективно увеличивает уровни высокочастотного шума.
hfvt = fvtool(df,[1 -1],1,'MagnitudeDisplay','zero-phase','Fs',Fs); legend(hfvt,'50th order FIR differentiator','Response of diff function');
Дифференцируйте использование различной функции. Добавьте нули, чтобы компенсировать недостающие выборки из-за различной операции.
v1 = diff(drift)/dt; a1 = diff(v1)/dt; v1 = [0; v1]; a1 = [0; 0; a1];
Дифференцируйтесь использование 50-го КИХ порядка фильтруют и компенсируют задержку.
D = mean(grpdelay(df)); % filter delay
v2 = filter(df,[drift; zeros(D,1)]);
v2 = v2(D+1:end);
a2 = filter(df,[v2; zeros(D,1)]);
a2 = a2(D+1:end);
v2 = v2/dt;
a2 = a2/dt^2;
Постройте несколько точек данных смещения пола. Постройте также несколько точек данных скорости и ускорения, как вычислено с diff и с 50-м КИХ-фильтром порядка. Заметьте, как шум был немного усилен в оценках скорости и в основном усилен в ускоряющих оценках, полученных с различным.
helperFilterIntroductionPlot2(t,drift,v1,v2,a1,a2)
Текучий фильтр интегратора является фильтром все-полюса с передаточной функцией![$H(Z) = 1/[1-cZ^{-1}]$](../../examples/signal/win64/FilterIntroductionExample_eq00578412725162239527.png)
, где
константа, которая должна быть меньшей, чем 1, чтобы гарантировать устойчивость фильтра. Не удивительно, что как подходы один, текучий интегратор приближается к инверсии различной передаточной функции. Примените текучий интегратор к ускорению и оценкам скорости, полученным в предыдущем разделе, чтобы возвратить скорость и дрейф соответственно. Используйте оценки, полученные с различной функцией, поскольку они являются более шумными.
Используйте текучий интегратор с
. Постройте ответ величины текучего фильтра интегратора. Заметьте, что фильтр действует как фильтр lowpass эффективно устранение высокочастотного шума.
close(hfvt)
fvtool(1,[1 -.999],'Fs',Fs)
Отфильтруйте скорость и ускорение с текучим интегратором.
v_original = v1;
a_original = a1;
d_leakyint = filter(1,[1 -0.999],v_original);
v_leakyint = filter(1,[1 -0.999],a_original);
% Multiply by time differential
d_leakyint = d_leakyint * dt;
v_leakyint = v_leakyint * dt;
Постройте смещение и оценки скорости и сравните с исходными сигналами v1 и a1.
helperFilterIntroductionPlot3(t,drift,d_leakyint,v_original,v_leakyint)
Можно также интегрировать сигнал с помощью функций cumtrapz и cumsum. Результаты будут похожи на полученных с текучим интегратором.
В этом примере вы узнали о линейных и нелинейных фильтрах фазы, и вы изучили, как компенсировать задержку фазы, введенную каждым типом фильтра. Вы также изучили, как применить фильтры, чтобы удалить нежелательные частотные составляющие из сигнала, и как проредить сигнал после ограничения его пропускной способности с фильтром lowpass. Наконец, вы изучили, как дифференцировать и интегрировать сигнал с помощью созданий цифровых фильтров. В примере вы также изучили, как использовать аналитические инструменты, чтобы посмотреть на ответ и групповую задержку ваших фильтров.
Для получения дополнительной информации о фильтре приложения видят Signal Processing Toolbox. Для получения дополнительной информации о том, как спроектировать цифровые фильтры, смотрите Практическое Введение в пример Создания цифровых фильтров.
Ссылки: Дж.Г. Проукис и Д. Г. Мэнолакис, "цифровая обработка сигналов. Принципы, алгоритмы и приложения", Prentice Hall, 1996.
С. Дж. Орфэнидис, "введение в обработку сигналов", Prentice Hall, 1996.
Следующие функции помощника используются в этом примере.