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

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

Спроектируйте 70-й порядок КИХ-фильтр lowpass с частотой среза 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 выполняет фильтрацию нулевой фазы путем обработки входных данных и в прямых и в обратных направлениях. Основной эффект состоит в том, что вы получаете искажение нулевой фазы, i.e., вы фильтруете данные с эквивалентным фильтром, который имеет постоянную задержку 0 выборок. Другие эффекты состоят в том, что вы получаете передаточную функцию фильтра, которая равняется величине в квадрате исходной передаточной функции фильтра и порядку фильтра, который удваивает порядок исходного фильтра.

Считайте сигнал ECG заданным в предыдущем разделе. Отфильтруйте этот сигнал с и незамедлительно компенсацию. Спроектируйте 7-й порядок БИХ lowpass эллиптический фильтр с частотой среза 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 отфильтрованный сигнал. Содержимое частоты выше 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\times 1400=2800$ Гц был бы достаточен, чтобы представлять сигнал правильно, вы однако, используют частоту дискретизации 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, чтобы достигнуть резкой метки частоты, маленькой неравномерности в полосе пропускания и относительно младшего разряда.

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

Дифференциация сигнала

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

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

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

Ссылки

Proakis, J. G. и Д. Г. Мэнолакис. Цифровая обработка сигналов: принципы, алгоритмы и приложения. Englewood Cliffs, NJ: Prentice Hall, 1996.

Orfanidis, S. J. Введение в обработку сигналов. Englewood Cliffs, NJ: Prentice Hall, 1996.

Приложение

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