Почему плохая идея фильтровать, обнуляя ячейки БПФ?

72

Очень легко отфильтровать сигнал, выполнив для него БПФ, обнулив некоторые ячейки, а затем выполнив IFFT. Например:

t = linspace(0, 1, 256, endpoint=False)
x = sin(2 * pi * 3 * t) + cos(2 * pi * 100 * t)
X = fft(x)
X[64:192] = 0
y = ifft(X)

Высокочастотный компонент полностью удаляется этим FFT-фильтром "кирпичной стены".

Но я слышал, что это не очень хороший метод для использования.

  • Почему это вообще плохая идея?
  • Есть ли обстоятельства, при которых это хороший или хороший выбор?

[ как предложено pichenettes ]

эндолиты
источник

Ответы:

74

Обнуление бинов в частотной области такое же, как умножение на прямоугольное окно в частотной области. Умножение на окно в частотной области аналогично циклической свертке путем преобразования этого окна во временной области. Преобразование прямоугольного окна - это функция Sinc ( ). Обратите внимание на то, что функция Sinc имеет много больших пульсаций и пульсаций, которые расширяют всю ширину апертуры временной области. Если фильтр во временной области, который может выводить все эти колебания (звонки), является «плохой идеей», то то же самое происходит и с нулевыми ячейками.sin(ωt)/ωt

Эти пульсации будут самыми большими для любого спектрального содержимого, которое является "между ячейками" или нецелым периодическим в ширине апертуры БПФ. Таким образом, если ваши исходные входные данные БПФ являются окном для каких-либо данных, которые в этом окне несколько непериодичны (например, большинство несинхронно дискретизированных сигналов «реального мира»), то эти конкретные артефакты будут создаваться с помощью нулевых бинов.

Другой способ взглянуть на это состоит в том, что каждый результирующий блок FFT представляет определенную частоту синусоидальной волны во временной области. Таким образом, обнуление ячейки даст тот же результат, что и вычитание этой синусоидальной волны, или, эквивалентно, добавление синусоидальной волны с точной центральной частотой ячейки БПФ, но с противоположной фазой. Обратите внимание, что если частота некоторого контента во временной области не является чисто целочисленной периодичностью по ширине БПФ, то попытка отменить ее, добавив инверсию ровно целочисленной периодической синусоидальной волны, приведет не к молчанию, а к чему-то более похожему нота "удара" (модулированная синусоидальная волна AM другой частоты). Опять же, вероятно, не то, что хотел.

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

hotpaw2
источник
3
У этого ответа есть хорошие вещи, но я бы предпочел больше сосредоточиться на эффекте Гиббса.
Джим Клей
4
Попытка получить ответ на эффект Гиббса уже задавалась здесь: dsp.stackexchange.com/questions/1144/…
hotpaw2
@ hotpaw2 Это хорошее объяснение. Тем не менее, мне нужна ссылка на это, и я нахожу некоторые трудности в ее идентификации. Это причина, по которой мы делаем фильтрацию во временной области, а не работаем в частотной области. (Кроме того, временная область может быть в реальном времени.) Однако никто, кажется, не начинает с того, что заявляет об этом!
Хью
Как это может быть связано с методом окна для дизайна фильтра?
Филип Пинто
Сравните преобразование окна фон Ханна (et.al.) с преобразованием любого прямоугольного окна. Гораздо лучше отклик фильтра в целом, особенно между бинами FFT в полосе останова. В общем, резкое обнуление бинов хуже, чем отсутствие обнуления вблизи переходов.
hotpaw2
3

Этот вопрос также давно смущал меня. Объяснение @ hotpaw2 хорошо. Возможно, вас заинтересует простой эксперимент с использованием matlab.

https://poweidsplearningpath.blogspot.com/2019/04/dftidft.html


обновленная информация.

Чтобы убедиться, что этот факт прост, нам просто нужно осторожно наблюдать спектр импульсной характеристики идеального (?) Полосового фильтра, который просто обнуляет ячейки БПФ. Почему мне нужно добавить наречие «осторожно»? Если мы просто используем тот же размер БПФ, чтобы наблюдать реакцию импульса, мы будем обмануты, как показано на рисунке 1 . Тем не менее, если мы добавим порядок ДПФ при наблюдении за выходным сигналом фильтра, то есть при заполнении нуля импульсной характеристикой, мы сможем найти так называемое явление Гиббса, рябь в частотной области, как показано на рисунке 2 .

На самом деле результат получается из оконного эффекта. Если вы хотите полностью понять проблему, обратитесь к главе 7.6 и главе 10.1-10.2 Библии о DSP (1). Подводя итог, можно отметить три ключевых момента.

  1. Размер окна и порядок DFT (FFT) полностью независимы. Не смешивайте их вместе.
  2. Свойства окна (тип / размер) доминируют над формой DTFT. (например, более широкие главные лепестки приводят к более широкой переходной полосе частотной характеристики.)
  3. DFT - это просто выборка DTFT в частотной области. Кроме того, чем выше порядок ДПФ, тем плотнее спектр ДПФ.

Итак, с помощью более плотного спектра на рис. 2 мы можем видеть сквозь маску идеального (поддельного) полосового фильтра.

введите описание изображения здесьОбманчиво Freq. Отклик.

введите описание изображения здесьФеномен Гиббса во Фреке. Отклик.

(1) Алан В. Оппенгейм и Рональд В. Шафер. 2009. Обработка сигналов с дискретным временем (3-е изд.). Прентис Холл Пресс, Аппер Седл Ривер, Нью-Джерси, США.

fps = 15;

LPF = 1;
HPF = 2;

n = -511:512;
n0 = 0;
imp = (n==n0);

NyquistF = 1/2*fps;

%% Ideal BPF
tmp_N = 512;
tmp_n = 0:1:tmp_N-1;
freq = ( n .* fps) ./ tmp_N;
F = fft(imp, tmp_N);  
F_bpf = IdealBandpassFilter(F, fps, LPF, HPF);
imp_rep =[real(ifft(F_bpf))'];

% Zero padding.
imp_rep2 =[zeros(1,2048) real(ifft(F_bpf))' zeros(1,2048)];

N = 2^nextpow2(length(imp_rep));
F = fft(imp_rep,N);
freq_step = fps/N;
freq = -fps/2:freq_step:fps/2-freq_step;
freq = freq(N/2+1:end)';

figure;
plot(freq,abs(F(1:N/2)));
xlabel('freq(Hz)');
ylabel('mag');
title('Mis leading Freq Response');


N = 2^nextpow2(length(imp_rep2));
F = fft(imp_rep2,N);
freq_step = fps/N;
freq = -fps/2:freq_step:fps/2-freq_step;
freq = freq(N/2+1:end)';

figure;
plot(freq,abs(F(1:N/2)));
xlabel('freq(Hz)');
ylabel('mag');
title('Zero Padding (DFT) with more points');

%% Function
function filered_signal = IdealBandpassFilter(input_signal, fs, w1, w2)

    N = length(input_signal);
    n = 0:1:N-1;
    freq = ( n .* fs) ./ N;

    filered_signal = zeros(N, 1);

    for i = 1:N
        if freq(i) > w1 & freq(i) < w2
            filered_signal(i) = input_signal(i);
        end

    end
end
По-вэй Хуан
источник
Можно ли это преобразовать в комментарий?
эндолит
Извините, мне не хватает репутации. Объяснение в URL тоже написано мной. Я просто хочу предоставить тестовый код, который может визуализировать плохое влияние пульсации.
По-Вэй Хуан
1

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

Обнуление бинов в FFT дает плохое разрешение после IFFT во временной области.

Итта Гутами
источник
однако есть вычислительные трудности для очень длинного сигнала, чтобы взять fft и затем ifft. Во избежание появления циттеров / звонков фильтрация сигнала должна проходить плавно от полосы пропускания к полосе останова.
Итта Гутами
«БПФ дает плохое временное разрешение» БПФ не дает временного разрешения, это преобразование спектральной области, и, как уже было сказано, дает только информацию о частотных составляющих сигнала.
EdParadox
Разрешение, предоставляемое БПФ, равно длине его окна. Все, что находится за пределами окна FFT, не разрешается как находящееся внутри окна FFT.
hotpaw2