Подтвердить что ты не робот

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

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

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

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

// 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) будут последним положительным частотный отсек.

Затем, чтобы уменьшить частоту bin 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 длины, и я всегда предоставляю все образцы, поэтому не требуется 0-дополнение.

// 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, где neg. значения частоты обнуляются перед синтезом.

Спасибо заранее.

4b9b3361

Ответ 1

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

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

Вот пример. Рассмотрим синусоидальный входной сигнал 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);

А вот график: IFFT of product

Сбой в начале блока - это совсем не то, чего мы ожидаем. Но если вы считаете fft(x), это имеет смысл. Дискретное преобразование Фурье предполагает, что сигнал является периодическим в блоке преобразования. Насколько известно ДПФ, мы просили преобразовать один период из этого: Aperiodic input to DFT

Это приводит к первому важному рассмотрению при фильтрации с помощью ДПФ: вы фактически реализуете круговую свертку, а не линейную свертку. Таким образом, "сбой" в первом графе на самом деле не является глюком, когда вы рассматриваете математику. Итак, возникает вопрос: существует ли способ работать с периодичностью? Ответ: дайте обработку перекрытия-сохранения. По сути, вы вычисляете продукты N-long, как указано выше, но сохраняете только 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: 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), вы получаете: hraw

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

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

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

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

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

Ответ 2

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

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

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

alt text

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

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

alt text

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

alt text

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

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

Ответ 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');

frequency response

Обратите внимание, что в моем коде я создаю пример значения 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'); 

frequency response

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

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

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

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

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

frequency response

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

Ответ 4

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

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