Фильтр нижних частот Matlab с использованием fft

Я реализовал простой фильтр нижних частот в Matlab, используя прямое и обратное преобразование Фурье. В принципе работает, но минимальное и максимальное значения отличаются от оригинала.

signal = data;
%% fourier spectrum
% number of elements in fft
NFFT = 1024;
% fft of data
Y = fft(signal,NFFT)/L;
% plot(freq_spectrum)

%% apply filter
fullw = zeros(1, numel(Y));
fullw( 1 : 20 ) = 1;
filteredData = Y.*fullw;

%% invers fft
iY = ifft(filteredData,NFFT);
% amplitude is in abs part
fY = abs(iY);
% use only the length of the original data
fY = fY(1:numel(signal));
filteredSignal = fY * NFFT; % correct maximum

clf; hold on;
plot(signal, 'g-')
plot(filteredSignal ,'b-')
hold off;

результирующее изображение выглядит так

введите здесь описание изображения

Что я делаю неправильно? Если я нормализую оба данных, отфильтрованный сигнал выглядит правильно.


matlab filtering fft lowpass-filter
person Matthias Pospiech    schedule 02.03.2015    source источник
comment
ваш фильтр должен быть симметричным, как и сигнал. почему вы ожидаете, что min и max не изменятся? нет причин.   -  person thang    schedule 02.03.2015
comment
Обратите внимание, что попытка применить фильтр кирпичной стены, подобный этому, в частотной области приведет к неприятным артефактам - вам нужно использовать гладкую функцию в частотной области (обычно функцию окна). Также обратите внимание, что усиление вашего фильтра не нормализовано, и, как отмечает @thang, ваш фильтр должен быть симметричным, иначе вы получите сложные выходные данные во временной области.   -  person Paul R    schedule 02.03.2015

Ответы (1)


arrow_upward
19
arrow_downward

Просто напомним себе, как MATLAB хранит частотный контент для Y = fft(y,N):

  • Y(1) - постоянное смещение
  • Y(2:N/2 + 1) - набор положительных частот
  • Y(N/2 + 2:end) — это набор отрицательных частот... (обычно мы наносим это слева от вертикальной оси)

Чтобы сделать настоящий фильтр нижних частот, мы должны сохранить как низкие положительные частоты , так и низкие отрицательные частоты.

Вот пример того, как это сделать с мультипликативным прямоугольным фильтром в частотной области, как вы сделали:

% make our noisy function
t = linspace(1,10,1024);
x = -(t-5).^2  + 2;
y = awgn(x,0.5); 
Y = fft(y,1024);

r = 20; % range of frequencies we want to preserve

rectangle = zeros(size(Y));
rectangle(1:r+1) = 1;               % preserve low +ve frequencies
y_half = ifft(Y.*rectangle,1024);   % +ve low-pass filtered signal
rectangle(end-r+1:end) = 1;         % preserve low -ve frequencies
y_rect = ifft(Y.*rectangle,1024);   % full low-pass filtered signal

hold on;
plot(t,y,'g--'); plot(t,x,'k','LineWidth',2); plot(t,y_half,'b','LineWidth',2); plot(t,y_rect,'r','LineWidth',2);
legend('noisy signal','true signal','+ve low-pass','full low-pass','Location','southwest')

введите здесь описание изображения

Полночастотный фитлер работает лучше, но вы заметите, что реконструкция немного «волнистая». Это связано с тем, что умножение с функцией прямоугольника в частотной области совпадает с сверткой с функцией sinc во времени. домен. Свертка с синусоидальной функцией заменяет каждую точку очень неравномерным средневзвешенным значением ее соседей, отсюда и эффект «волны».

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

gauss = zeros(size(Y));
sigma = 8;                           % just a guess for a range of ~20
gauss(1:r+1) = exp(-(1:r+1).^ 2 / (2 * sigma ^ 2));  % +ve frequencies
gauss(end-r+1:end) = fliplr(gauss(2:r+1));           % -ve frequencies
y_gauss = ifft(Y.*gauss,1024);

hold on;
plot(t,x,'k','LineWidth',2); plot(t,y_rect,'r','LineWidth',2); plot(t,y_gauss,'c','LineWidth',2);
legend('true signal','full low-pass','gaussian','Location','southwest')

введите здесь описание изображения

Как видите, реконструкция получается намного лучше.

person eigenchris    schedule 02.03.2015
comment
Очень красивое объяснение. Признаюсь, я склонен забывать, что пространство Фурье свернуто и фильтрацию нужно делать с обеих сторон. - person Matthias Pospiech; 03.03.2015