Компенсация задержки, введенной фильтрацией

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

Задержка фильтра, которая является постоянной на всех частотах, может быть легко скомпенсирована сдвигом сигнала во времени. Фильтры FIR обычно имеют постоянную задержку. С другой стороны, задержка, которая изменяется с частотой, вызывает фазовые искажения и может значительно изменить форму сигнала. Компенсация частотно-зависимой задержки не так тривиальна, как для случая постоянной задержки. БИХ-фильтры вводят частотно-зависимую задержку.

Компенсация постоянной задержки фильтра

Как упоминалось выше, можно измерить группу задержки фильтра, чтобы убедиться, что он является постоянной функцией частоты. Вы можете использовать grpdelay функция для измерения задержки D фильтра и компенсации этой задержки путем добавления D нулей к входному сигналу и сдвига выходного сигнала во времени на D выборок.

Рассмотрите шумной сигнал электрокардиограммы, который вы хотите отфильтровать, чтобы удалить высокочастотный шум выше 75 Гц. Вы хотите применить фильтр нижних частот FIR и компенсировать задержку фильтра, чтобы шумные и отфильтрованные сигналы были правильно выровнены и могли быть нанесены друг на друга для сравнения.

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

Сконструировать фильтр FIR нижних частот 70-го порядка с частотой отсечки 75 Гц.

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

plot(t,x,t,y,'linewidth',1.5)
title('Filtered Waveforms')
xlabel('Time (s)')
legend('Original Noisy Signal','Filtered Signal')
grid on
axis tight

Компенсация частотно-зависимой задержки

Частотно-зависимая задержка вызывает фазовые искажения в сигнале. Компенсация этого типа задержки не так тривиальна, как для случая постоянной задержки. Если приложение поддерживает автономную обработку, можно удалить частотно-зависимую задержку, реализовав нулевую фазовую фильтрацию с помощью filtfilt функция. filtfilt выполняет нулевую фазовую фильтрацию путем обработки входных данных как в прямом, так и в обратном направлениях. Основной эффект состоит в том, что получается искажение нулевой фазы, т.е. данные фильтруются эквивалентным фильтром с постоянной задержкой 0 выборок. Другие эффекты заключаются в том, что получается передаточная функция фильтра, которая равна квадратной величине исходной передаточной функции фильтра, и порядок фильтров, который вдвое превышает порядок исходного фильтра.

Рассмотрим сигнал ЭКГ, определенный в предыдущем разделе. Фильтрация этого сигнала с компенсацией задержки и без нее. Проектирование низкочастотного БИХ-эллиптического фильтра 7-го порядка с частотой отсечки 75 Гц.

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);    % Nonlinear phase filter - no delay compensation
y2 = filtfilt(df,x);  % Zero-phase implementation - delay compensation

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')
xlim([0.25 0.55])
grid on

Фильтрация по нулевой фазе - это отличный инструмент, если приложение позволяет выполнять некаузальные операции прямой/обратной фильтрации, а также изменять отклик фильтра на квадрат исходного отклика.

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

Удаление нежелательного спектрального содержимого из сигнала

Фильтры обычно используются для удаления нежелательного спектрального содержимого из сигнала. Для этого можно выбрать из множества фильтров. Фильтр нижних частот выбирается при удалении высокочастотного содержимого или фильтр верхних частот при удалении низкочастотного содержимого. Также можно выбрать полосовой фильтр для удаления низкочастотного и высокочастотного содержимого, оставляя интактным промежуточный диапазон частот. Вы выбираете полосовой фильтр, когда хотите удалить частоты в данной полосе.

Рассмотрим аудиосигнал с гулом линии питания и белым шумом. Гул линии питания вызван тональным сигналом 60 Гц. Белый шум - это сигнал, который существует по всей полосе пропускания звука.

Загрузите звуковой сигнал. Укажите частоту дискретизации 44,1 кГц.

Fs = 44100;
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'})

Можно сначала удалить как можно больше спектрального содержимого белого шума с помощью фильтра нижних частот. Полоса пропускания фильтра должна быть установлена на значение, которое обеспечивает хороший компромисс между снижением шума и ухудшением качества звука из-за потери высокочастотного содержимого. Применение фильтра нижних частот перед удалением hum 60 Гц очень удобно, так как вы сможете уменьшить диапазон сигнала. Сигнал более низкой скорости позволит сконструировать более резкий и узкий 60 Гц полосовой фильтр с меньшим порядком фильтров.

Спроектировать фильтр нижних частот с частотой полосы пропускания 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

df = designfilt('lowpassfir','PassbandFrequency',Fp,...
                'StopbandFrequency',Fst,'PassbandRipple',Ap,...
                'StopbandAttenuation',Ast,'SampleRate',Fs);

Проанализируйте отклик фильтра.

fvtool(df,'Fs',Fs,'FrequencyScale','log',...
  'FrequencyRange','Specify freq. vector','FrequencyVector',F)

Фильтрация данных и компенсация задержки.

D = mean(grpdelay(df)); % Filter delay
ylp = filter(df,[y; zeros(D,1)]);
ylp = ylp(D+1:end);

Посмотрите на спектр отфильтрованного сигнала нижних частот. Удалено частотное содержание выше 1400 Гц.

[Plp,Flp] = pwelch(ylp,ones(8192,1),8192/2,8192,Fs,'power');
helperFilterIntroductionPlot1(F,P,Flp,Plp,...
  {'Original signal','Lowpass filtered signal'})

На приведенном выше графике спектра мощности видно, что максимальное не пренебрежимо малое частотное содержание фильтруемого сигнала нижних частот составляет 1400 Гц. По теореме дискретизации, частота дискретизации $2\mathrm{\times }\mathrm{1400\; =}$ 2800 Гц было бы достаточно, чтобы представить сигнал правильно, вы, однако, используете частоту дискретизации 44100 Гц, которая является пустой, так как вам нужно будет обработать больше образцов, чем необходимо. Можно уменьшить частоту дискретизации и вычислительную нагрузку, уменьшив количество выборок, которые необходимо обработать. Более низкая частота дискретизации также позволит сконструировать более резкий и узкий полосовой фильтр, необходимый для удаления шума 60 Гц, с меньшим порядком фильтров.

Для получения частоты дискретизации, равной Fs/10 = 4,41 кГц, выполните обратную выборку отфильтрованного сигнала нижних частот с коэффициентом 10. Постройте график спектра сигнала до и после понижающей дискретизации.

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 Гц с помощью фильтра bandstop IIR. Пусть полоса останова имеет ширину 4 Гц с центром 60 Гц. Выберите фильтр IIR, чтобы получить резкий скачок частоты, небольшую пульсацию полосы пропускания и относительно низкий порядок.

df = designfilt('bandstopiir','PassbandFrequency1',55,...
               'StopbandFrequency1',58,'StopbandFrequency2',62,...
               'PassbandFrequency2',65,'PassbandRipple1',1,...
               'StopbandAttenuation',60,'PassbandRipple2',1,...
               'SampleRate',Fs,'DesignMethod','ellip');                          

Проанализируйте отклик величины.

fvtool(df,'Fs',Fs,'FrequencyScale','log',...
  'FrequencyRange','Specify freq. vector','FrequencyVector',Fds(Fds>F(2)))

Во избежание фазовых искажений выполните фильтрацию по нулевой фазе. Используйте filtfilt для обработки данных.

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');
helperFilterIntroductionPlot1(F,P,Ffinal,Pfinal,...
  {'Original signal','Final filtered signal'})

Если на компьютере установлена звуковая карта, можно прослушать сигнал до и после обработки. Как упоминалось выше, конечный результат состоит в том, что вы эффективно ослабили гул 60 Гц и высокочастотный шум в аудиофайле.

% To play the original signal, uncomment the next two lines
% hplayer = audioplayer(y, Fs);
% play(hplayer)

% To play the noise-reduced signal, uncomment the next two lines
% hplayer = audioplayer(yf, Fs);
% play(hplayer);

Дифференцирование сигнала

MATLAB diff функция дифференцирует сигнал с тем недостатком, что можно потенциально увеличить уровни шума на выходе. Лучшим вариантом является использование фильтра дифференциатора, который действует как дифференциатор в интересующей полосе и как аттенюатор на всех других частотах, эффективно устраняя высокочастотный шум.

В качестве примера проанализируйте скорость смещения перекрытия здания во время землетрясения. Измерения смещения или дрейфа регистрировались на первом этаже трехэтажной испытательной конструкции в условиях землетрясения и сохранялись в файле quakedrift.mat. Длина вектора данных равна 10e3, частота дискретизации равна 1 кГц, а единицы измерения составляют см..

Дифференцировать данные смещения для получения оценок скорости и ускорения перекрытия здания во время землетрясения. Сравните результаты с помощью фильтра дифференциатора diff и FIR.

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

diff функция может рассматриваться как фильтр FIR первого порядка с ответом $$H(Z)\mathrm{}={}^{\mathrm{}}$$1-Z-1. Используйте FVTool, чтобы сравнить отклик величины дифференцирующего FIR фильтра 50-го порядка и отклик diff функция. Очевидно, что оба ответа эквивалентны в области полосы пропускания (от 0 до 100 Гц). Однако в области стоп-полосы фильтр 50-го порядка ослабляет компоненты, в то время как отклик diff усиливает компоненты. Это эффективно увеличивает уровни высокочастотного шума.

hfvt = fvtool(df,[1 -1],1,'MagnitudeDisplay','zero-phase','Fs',Fs);
legend(hfvt,'50th order FIR differentiator','Response of diff function');

Дифференцировать с помощью diff функция. Добавьте нули для компенсации недостающих выборок из-за операции разбиения.

v1 = diff(drift)/dt;
a1 = diff(v1)/dt;

v1 = [0; v1];
a1 = [0; 0; a1];

Дифференцировать с помощью фильтра FIR 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;

Постройте график нескольких точек данных смещения пола. Постройте график также нескольких точек данных скорости и ускорения, вычисленных с помощью фильтра FIR 50-го порядка. Обратите внимание, как шум был слегка усилен в оценках скорости и в значительной степени усилен в оценках ускорения, полученных с помощью diff.

helperFilterIntroductionPlot2(t,drift,v1,v2,a1,a2)

Интеграция сигнала

Фильтр интегратора с утечкой представляет собой общеполюсный фильтр с передаточной функцией $$H(Z)\mathrm{}=1/[{}^{\mathrm{}}$$1-cZ-1], $$$$где c - постоянная, которая должна быть меньше 1 для обеспечения стабильности фильтра. Неудивительно, что $$$$когда c приближается к единице, протекающий интегратор приближается к обратному diff передаточная функция. Примените интегратор утечки к оценкам ускорения и скорости, полученным в предыдущем разделе, чтобы вернуть скорость и дрейф соответственно. Использовать оценки, полученные с помощью the diff функция, так как они более шумные.

Используйте протекающий интегратор с $$a\mathrm{=}\mathrm{}\mathrm{}\mathrm{}$$0,999. Постройте график амплитудной характеристики фильтра интегратора с утечкой. Фильтр действует как фильтр нижних частот, эффективно устраняя высокочастотный шум.

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

d_leakyint = d_leakyint * dt;
v_leakyint = v_leakyint * dt;

Постройте график оценок смещения и скорости и сравните их с исходными сигналами.

helperFilterIntroductionPlot3(t,drift,d_leakyint,v_original,v_leakyint)

Также можно интегрировать сигнал с помощью cumsum и cumtrapz функции. Результаты будут аналогичны результатам, полученным с помощью протекающего интегратора.

Заключения

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

Для получения дополнительной информации о приложениях фильтров см. документацию по Toolbox™ обработки сигналов. Дополнительные сведения о проектировании цифровых фильтров см. в разделе Практическое введение в проектирование цифровых фильтров.

Ссылки

Проакис, Дж. Г. и Д. Г. Манолакис. Цифровая обработка сигналов: принципы, алгоритмы и приложения. Энглвуд Клиффс, Нью-Джерси: Прентис-Холл, 1996.

Орфанидис, С. Дж. Введение в обработку сигналов. Энглвуд Клиффс, Нью-Джерси: Прентис-Холл, 1996.

Приложение

В этом примере используются следующие вспомогательные функции: