DSP - фильтрация в частотной области через FFT

Question

Я немного поигрался с реализацией FFT Exocortex, но у меня возникли некоторые проблемы.

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

Позвольте мне привести пример выходного буфера FFT с 4 выборками:

// Bin 0 (DC)
FFTOut[0] = 0.0000610351563
FFTOut[1] = 0.0

// Bin 1
FFTOut[2] = 0.000331878662
FFTOut[3] = 0.000629425049

// Bin 2
FFTOut[4] = -0.0000381469727
FFTOut[5] =  0.0

// Bin 3, this is the first and only negative frequency bin.
FFTOut[6] =  0.000331878662
FFTOut[7] = -0.000629425049

Выходные данные состоят из пар чисел с плавающей точкой, каждая из которых представляет действительную и воображаемую части одного бина. Таким образом, bin 0 (индексы массива 0, 1) будет представлять действительную и мнимую части частоты постоянного тока. Как вы можете видеть, бины 1 и 3 имеют одинаковые значения (за исключением знака части Im), поэтому я предполагаю, что бин 3 является первой отрицательной частотой, и, наконец, индексы (4, 5) будут последними положительными частотный блок.

Затем, чтобы ослабить частотный интервал 1, я делаю это:

// Attenuate the 'positive' bin
FFTOut[2] *= 0.5;
FFTOut[3] *= 0.5;

// Attenuate its corresponding negative bin.
FFTOut[6] *= 0.5;
FFTOut[7] *= 0.5;

Для реальных тестов я использую БПФ длиной 1024 и всегда предоставляю все образцы, поэтому заполнение нулями не требуется.

// Attenuate
var halfSize = fftWindowLength / 2;
float leftFreq = 0f;
float rightFreq = 22050f; 
for( var c = 1; c < halfSize; c++ )
{
    var freq = c * (44100d / halfSize);

    // Calc. positive and negative frequency indexes.
    var k = c * 2;
    var nk = (fftWindowLength - c) * 2;

    // This kind of attenuation corresponds to a high-pass filter.
    // The attenuation at the transition band is linearly applied, could
    // this be the cause of the distortion of low frequencies?
    var attn = (freq < leftFreq) ? 
                    0 : 
                    (freq < rightFreq) ? 
                        ((freq - leftFreq) / (rightFreq - leftFreq)) :
                        1;

    // Attenuate positive and negative bins.
    mFFTOut[ k ] *= (float)attn;
    mFFTOut[ k + 1 ] *= (float)attn;
    mFFTOut[ nk ] *= (float)attn;
    mFFTOut[ nk + 1 ] *= (float)attn;
}

Очевидно, что я делаю что-то не так, но не могу понять, что.

Я не хочу использовать выход FFT в качестве средства для генерации набора коэффициентов FIR, поскольку я пытаюсь реализовать очень простой динамический эквалайзер.

Какой правильный способ фильтрации в частотной области? чего мне не хватает?

Кроме того, действительно ли необходимо также ослаблять отрицательные частоты? Я видел реализацию FFT, где нег. Значения частоты обнуляются перед синтезом.

Заранее спасибо.

Answer 1

Есть две проблемы: способ использования БПФ и конкретный фильтр.

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

Вот пример. Рассмотрим синусоидальный входной сигнал x с 1,5 циклами в периоде и простой фильтр нижних частот h. В синтаксисе Matlab / Octave:

N = 1024;
n = (1:N)'-1; %'# define the time index
x = sin(2*pi*1.5*n/N); %# input with 1.5 cycles per 1024 points
h = hanning(129) .* sinc(0.25*(-64:1:64)'); %'# windowed sinc LPF, Fc = pi/4
h = [h./sum(h)]; %# normalize DC gain

y = ifft(fft(x) .* fft(h,N)); %# inverse FT of product of sampled spectra
y = real(y); %# due to numerical error, y has a tiny imaginary part
%# Depending on your FT/IFT implementation, might have to scale by N or 1/N here
plot(y);

А вот график:

Ошибка в начале блока совсем не та, которую мы ожидаем. Но если учесть fft(x), это имеет смысл. Дискретное преобразование Фурье предполагает, что сигнал является периодическим в пределах блока преобразования. Насколько известно DFT, мы попросили преобразовать один период этого:

Это приводит к первому важному соображению при фильтрации с помощью ДПФ: вы фактически реализуете круговую свертку , а не линейную свертку. Таким образом, «сбой» на первом графике на самом деле не сбой, если учесть математику. Тогда возникает вопрос: есть ли способ обойти периодичность? Ответ - да: используйте обработка с сохранением перекрытия . По сути, вы рассчитываете N-длинные продукты, как указано выше, но сохраняете только N / 2 балла.

Nproc = 512;
xproc = zeros(2*Nproc,1); %# initialize temp buffer
idx = 1:Nproc; %# initialize half-buffer index
ycorrect = zeros(2*Nproc,1); %# initialize destination
for ctr = 1:(length(x)/Nproc) %# iterate over x 512 points at a time
    xproc(1:Nproc) = xproc((Nproc+1):end); %# shift 2nd half of last iteration to 1st half of this iteration
    xproc((Nproc+1):end) = x(idx); %# fill 2nd half of this iteration with new data
    yproc = ifft(fft(xproc) .* fft(h,2*Nproc)); %# calculate new buffer
    ycorrect(idx) = real(yproc((Nproc+1):end)); %# keep 2nd half of new buffer
    idx = idx + Nproc; %# step half-buffer index
end

А вот график ycorrect:

Эта картина имеет смысл - мы ожидаем переходный процесс запуска от фильтра, а затем результат устанавливается в синусоидальный отклик установившегося состояния. Обратите внимание, что теперь x может быть произвольно длинным. Ограничение Nproc > 2*min(length(x),length(h)).

Теперь перейдем ко второму вопросу: конкретному фильтру. В вашем цикле вы создаете фильтр, спектр которого по существу равен H = [0 (1:511)/512 1 (511:-1:1)/512]';. Если вы делаете hraw = real(ifft(H)); plot(hraw), вы получаете:

Трудно увидеть, но есть куча ненулевых точек на дальнем левом краю графика, а затем еще кучка на дальнем правом краю. Используя встроенную в Octave функцию freqz, чтобы посмотреть на частотную характеристику, которую мы видим (выполнив freqz(hraw)):

Отклик величины имеет много пульсаций от огибающей верхних частот до нуля. Опять же, периодичность, присущая DFT, работает. Что касается DFT, hraw повторяется снова и снова. Но если взять один период hraw, как freqz, его спектр будет сильно отличаться от периодического варианта.

Итак, давайте определим новый сигнал: hrot = [hraw(513:end) ; hraw(1:512)]; Мы просто вращаем необработанный вывод DFT, чтобы сделать его непрерывным внутри блока. Теперь давайте посмотрим на частотную характеристику, используя freqz(hrot):

Намного лучше. Желаемый конверт там, без всякой ряби. Конечно, реализация теперь не так проста, вам нужно сделать полное комплексное умножение на fft(hrot) вместо простого масштабирования каждого сложного бина, но по крайней мере вы получите правильный ответ.

Обратите внимание, что для скорости вы обычно предварительно вычисляете DFT для дополненного h, я оставил его в цикле, чтобы было проще сравнивать его с оригиналом.

Answer 2

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

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

Например, взяв несколько разумных чисел: я начинаю с формы волны с частотой 27,5 Гц (A0 на фортепиано), оцифрованной с частотой 44100 Гц, это будет выглядеть так (где красная часть имеет длину 1024 семпла):

альтернативный текст http://i48.tinypic.com/zim802.png

Итак, сначала мы начнем с низких частот 40 Гц. Таким образом, поскольку исходная частота меньше 40 Гц, фильтр нижних частот с частотой среза 40 Гц не должен иметь никакого эффекта, и мы получим выход, который почти точно соответствует входу. Правильно? Неправильно, неправильно, неправильно - и это в основном суть вашей проблемы. Проблема в том, что для коротких участков идея с частотой 27,5 Гц четко не определена и не может быть хорошо представлена ​​в DFT.

То, что 27,5 Гц не имеют особого значения в коротком сегменте, можно увидеть, посмотрев на ДПФ на рисунке ниже. Обратите внимание, что хотя DFT более длинного сегмента (черные точки) показывает пик на 27,5 Гц, короткий (красные точки) не имеет.

альтернативный текст http://i50.tinypic.com/14w6luw.png

Очевидно, что при фильтрации ниже 40 Гц будет просто фиксироваться смещение постоянного тока, а результат фильтра нижних частот 40 Гц показан зеленым цветом ниже.

альтернативный текст http://i48.tinypic.com/2vao21w.png

Синяя кривая (взятая с отсечкой 200 Гц) начинает намного лучше совпадать. Но обратите внимание, что это не низкие частоты, которые делают его хорошо совпадают, а включение высоких частот. Только когда мы включим каждую возможную частоту в короткий сегмент, вплоть до 22 кГц, мы наконец получим хорошее представление об исходной синусоиде.

Причина всего этого заключается в том, что небольшой сегмент синусоидальной волны 27,5 Гц равен , а не синусоидальной частоте 27,5 Гц, и его DFT не имеет ничего общего с 27,5 Гц.

Answer 3

Вы ослабляете значение выборки частоты постоянного тока до нуля? Похоже, что вы не ослабляете это вообще в своем примере. Поскольку вы реализуете фильтр верхних частот, вам также необходимо установить значение DC на ноль.

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

Вот пример в MATLAB / Octave для демонстрации того, что может происходить:

N = 32;
os = 4;
Fs = 1000;
X = [ones(1,4) linspace(1,0,8) zeros(1,3) 1 zeros(1,4) linspace(0,1,8) ones(1,4)];
x = ifftshift(ifft(X));
Xos = fft(x, N*os);
f1 = linspace(-Fs/2, Fs/2-Fs/N, N);
f2 = linspace(-Fs/2, Fs/2-Fs/(N*os), N*os);

hold off;
plot(f2, abs(Xos), '-o');
hold on;
grid on;
plot(f1, abs(X), '-ro');
hold off;
xlabel('Frequency (Hz)');
ylabel('Magnitude');

АЧХ http://www.freeimagehosting.net/uploads/e10109e535.png

Обратите внимание, что в моем коде я создаю пример значения DC, отличного от нуля, с последующим резким изменением на ноль и затем увеличением. Затем я беру IFFT для преобразования во временную область. Затем я выполняю заполненное нулями fft (которое выполняется автоматически MATLAB, когда вы передаете размер fft больше входного сигнала) для этого сигнала во временной области. Заполнение нулями во временной области приводит к интерполяции в частотной области. Используя это, мы можем видеть, как фильтр будет реагировать между выборками фильтра.

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

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

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

N = 32;
os = 4;
Fs = 1000;
X = [ones(1,8) zeros(1,16) ones(1,8)];
x = ifftshift(ifft(X));
Xos = fft(x, N*os);
f1 = linspace(-Fs/2, Fs/2-Fs/N, N);
f2 = linspace(-Fs/2, Fs/2-Fs/(N*os), N*os);

hold off;
plot(f2, abs(Xos), '-o');
hold on;
grid on;
plot(f1, abs(X), '-ro');
hold off;
xlabel('Frequency (Hz)');
ylabel('Magnitude'); 

АЧХ http://www.freeimagehosting.net/uploads/c86f5f1700.png

Обратите внимание, что существует много колебаний, вызванных резкими изменениями.

БПФ или дискретное преобразование Фурье является выборочной версией преобразования Фурье. Преобразование Фурье применяется к сигналу в непрерывном диапазоне от бесконечности до бесконечности, тогда как ДПФ применяется к конечному числу выборок. Это в действительности приводит к квадратному оконному отображению (усечению) во временной области при использовании ДПФ, поскольку мы имеем дело только с конечным числом выборок. К сожалению, ДПФ прямоугольной волны является функцией типа sinc (sin (x) / x).

Проблема с резкими переходами в вашем фильтре (быстрый переход от 0 до 1 в одном примере) состоит в том, что он имеет очень длинный отклик во временной области, который усекается квадратным окном. Таким образом, чтобы минимизировать эту проблему, мы можем умножить сигнал во временной области на более постепенное окно. Если мы умножим окно Ханнинга, добавив строку:

x = x .* hanning(1,N).';

после принятия IFFT мы получаем такой ответ:

АЧХ http://www.freeimagehosting.net/uploads/944da9dd93.png

Поэтому я бы порекомендовал попытаться реализовать метод проектирования окон, поскольку он довольно прост (есть лучшие способы, но они усложняются). Поскольку вы внедряете эквалайзер, я предполагаю, что вы хотите иметь возможность изменять затухания на лету, поэтому я бы посоветовал рассчитывать и хранить фильтр в частотной области при каждом изменении параметров, а затем вы можете просто применить его к каждому входному аудиобуферу, взяв fft входного буфера, умножив его на выборки фильтра в частотной области, а затем выполнив ifft, чтобы вернуться во временную область. Это будет намного эффективнее, чем все ветвления, которые вы делаете для каждого образца.

Answer 4

Во-первых, о нормализации: это известная (не) проблема. DFT / IDFT потребует коэффициент 1 / sqrt (N) (кроме стандартных косинус / синусоидальных факторов) в каждом из них (прямой обратный), чтобы сделать их симметричными и действительно обратимыми. Другая возможность состоит в том, чтобы разделить одно из них (прямое или обратное) на N , это вопрос удобства и вкуса. Часто подпрограммы FFT не выполняют эту нормализацию, пользователь должен знать об этом и нормализовать, как он предпочитает. См

Второе: в (скажем) 16-точечном ДПФ то, что вы называете bin 0 , будет соответствовать нулевой частоте (DC), bin 1 low freq ... bin 4 средняя частота, bin 8 до самой высокой частоты и bin 9 ... 15 до "отрицательных частот". В нашем примере bin 1 на самом деле является как низкой, так и средней частотой. Помимо этого соображения, нет ничего концептуально неправильного в вашем "уравнении". Я не понимаю, что вы подразумеваете под «сигнал искажается на низких частотах» . Как вы это наблюдаете?