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

Цифровые фильтры вводят задержку в вашем сигнале. В зависимости от характеристик фильтра, задержка может быть постоянной по всем частотам, или она может варьироваться с частотой. Тип задержки определяет действия, которые вы должны предпринять, чтобы компенсировать его. The 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

Проектируйте lowpass конечная импульсная характеристика 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 выборки. Другие эффекты заключаются в том, что вы получаете передаточную функцию фильтра, которая равна квадрату величины исходной передаточной функции фильтра, и порядок фильтра, который вдвое превышает порядок исходного фильтра.

Учитывайте сигнал ЭКГ, заданный в предыдущем разделе. Пропустите этот сигнал с компенсацией и без нее. Проектируйте lowpass БИХ фильтр 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

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

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

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

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

Рассмотрим аудиосигнал, который имеет гул линии питания и белый шум. Гул линии питания вызван тоном 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'})

Можно сначала удалить как можно больше спектрального содержимое белого шума с помощью 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

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

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

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

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

Понижайте частоту фильтрованного сигнала lowpass в 10 раз, чтобы получить частоту дискретизации Fs/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 Гц. Выберите БИХ фильтр, чтобы получить резкую частотную надрез, маленькие неравномерности в полосе пропускания и относительно низкий порядок.

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 и конечной импульсной характеристики дифференциатора.

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

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

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

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

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-го порядка. Заметьте, как шум был слегка усилен в оценках скорости и в значительной степени усилен в оценках ускорения, полученных с помощью diff.

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

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

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

Используйте утечку интегратора с $$a=0.999$$. Постройте график амплитудной характеристики фильтра утечки интегратора. Фильтр действует как lowpass, эффективно устраняя высокочастотный шум.

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 функций. Результаты будут аналогичны результатам, полученным с утечкой интегратора.

Заключения

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

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

Ссылки

Проакис, Дж. Г., и Д. Г. Манолакис. Цифровая обработка сигналов: принципы, алгоритмы и приложения. Englewood Cliffs, Нью-Джерси: Prentice Hall, 1996.

Орфанидис, С. Дж. Введение в обработку сигналов. Englewood Cliffs, Нью-Джерси: Prentice Hall, 1996.

Приложение

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