DSP — Фильтрация в частотной области посредством БПФ

Question

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

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

Позвольте мне привести пример выходного буфера БПФ с 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

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

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

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

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

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

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

Answer 1

Есть два вопроса: то, как вы используете FFT и конкретный фильтр.

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

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

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, мы просили преобразование одного периода этого:

Это приводит к первому важным рассмотрении при фильтрации с DFTS: вы на самом деле реализуете Круговая свертка, не линейная свертка. Таким образом, «глюк» на первом графике не совпадает с глюдом, когда вы считаете математику. Так что тогда вопрос становится: есть ли способ работать вокруг периодичности? Ответ да: использовать Перекрытие обработки. Отказ По сути, вы рассчитываете N-Dange Products, как указано выше, но только держите 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), ты получаешь:

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

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

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

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

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

Answer 2

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

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

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

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

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

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

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

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

Answer 3

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

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

Вот пример в 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');

Обратите внимание, что в моем коде я создаю пример значения постоянного тока, ненулена, а затем резкое изменение на ноль, а затем поднимается. Затем я возьму преобразование IFFT, чтобы превратиться в момент времени. Затем я выполняю нулевой FFT 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'); 

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

Преобразование FFT или дискретных Фурье - это выборка преобразования Фурье. Преобразование Фурье применяется к сигналу по поводу непрерывного диапазона до бесконечности, в то время как DFT наносится над конечным количеством образцов. Это в результате чего приводит к квадратному оженному (усечению) в момент времени при использовании DFT, поскольку мы имеем дело только с конечным количеством образцов. К сожалению, DFT квадратной волны - функция типа SINC (SIN (X) / X).

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

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

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

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

Answer 4

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

Во-вторых: в A (говорят) 16 пунктов DFT, что вы называете Bin 0. будет соответствовать нулевой частоте (DC), Bin 1. Низкий Freq ... Bin 4 Средний Freq, Bin 8. на самую высокую частоту и BINS 9 ... 15 на «отрицательные частоты». В тебе пример, то, то, Bin 1. на самом деле как низкая частота, так и средняя частота. Помимо этого рассмотрения, в вашем «выравнивании» нет ничего особенного. Я не понимаю, что вы подразумеваете под «Сигнал искажен на низких частотах». Отказ Как вы считаете этим?