DSP - Просачивание частотной области через FFT
Я играл вокруг немного с реализацией Экзокоры FFT, но у меня есть некоторые проблемы.
Каждый раз, когда я изменяю амплитуды мусорных ведер частоты прежде, чем назвать 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
Вывод состоит из пар плаваний, каждый представляющий реальные и мнимые части единственного мусорного ведра. Так, мусорное ведро 0 (индексы массива 0, 1) представило бы действительные и мнимые части частоты DC. Как Вы видите, мусорные ведра 1, и 3 у обоих есть те же значения, (за исключением знака меня часть), таким образом, я предполагаю, что мусорное ведро 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;
Для фактических тестов я использую FFT с 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, где отрицательные значения частоты обнуляются перед синтезом.
Заранее спасибо.
4 ответа
Есть две проблемы: способ использования БПФ и конкретный фильтр.
Фильтрация традиционно реализуется как свертка во временной области. Вы правы, что перемножение спектров входного сигнала и сигнала фильтра эквивалентно. Однако, когда вы используете дискретное преобразование Фурье (ДПФ) (реализованное с помощью алгоритма быстрого преобразования Фурье для скорости), вы фактически вычисляете выборочную версию истинного спектра. Это имеет множество последствий, но наиболее важным для фильтрации является то, что сигнал во временной области является периодическим.
Вот пример. Рассмотрим синусоидальный входной сигнал 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), то это имеет смысл. Дискретное преобразование Фурье предполагает, что сигнал периодичен в пределах блока преобразования. Насколько известно ДПФ, мы попросили преобразовать один период этого сигнала:
Это приводит к первому важному соображению при фильтрации с помощью ДПФ: на самом деле вы реализуете круговую свертку, а не линейную. Поэтому "глюк" на первом графике на самом деле не является глюком, если рассматривать математику. Тогда возникает вопрос: есть ли способ обойти периодичность? Ответ - да: используйте обработку overlap-save. По сути, вы вычисляете 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)):
Частотная характеристика имеет много пульсаций от огибающей высоких частот до нуля. Снова работает периодичность, присущая ДПФ. Что касается ДПФ, hraw повторяется снова и снова. Но если взять один период hraw, как это делает freqz, то его спектр сильно отличается от спектра периодической версии.
Итак, определим новый сигнал: hrot = [hraw(513:end) ; hraw(1:512)]; Мы просто поворачиваем необработанный выход DFT, чтобы сделать его непрерывным внутри блока. Теперь давайте посмотрим на частотную характеристику, используя freqz(hrot):
Намного лучше. Желаемая огибающая есть, но без всех этих пульсаций. Конечно, теперь реализация не так проста, вам придется делать полное комплексное умножение на fft(hrot), а не просто масштабировать каждый комплексный бином, но, по крайней мере, вы получите правильный ответ.
Обратите внимание, что для скорости вы обычно предварительно вычисляете ДПФ заполненного h, я оставил его в цикле, чтобы легче было сравнить с оригиналом.
Во-первых, о нормализации: это известная (не) проблема. Для DFT/IDFT потребуется фактор 1/sqrt(N) (помимо стандартных факторов косинус/синус) в каждом из них (прямом и обратном), чтобы сделать их симметричными и действительно инвертируемыми. Другой возможностью является деление одного из них (прямого или обратного) на N, это вопрос удобства и вкуса. Часто процедуры БПФ не выполняют эту нормализацию, ожидается, что пользователь будет знать об этом и выполнять нормализацию по своему усмотрению. См.
Второе: в (скажем) 16-точечном ДПФ то, что вы называете bin 0, будет соответствовать нулевой частоте (DC), bin 1 низкой частоте... bin 4 средней частоте, bin 8 - самой высокой частоте и bin 9...15 - "отрицательным частотам". Таким образом, в вашем примере бин 1 фактически является и низкой, и средней частотой. Кроме этого соображения, в вашей "эквализации" нет ничего концептуально неправильного. Я не понимаю, что вы имеете в виду под "сигнал искажается на низких частотах". Как вы это наблюдаете?
Вы уменьшаете значение выборки частоты постоянного тока до нуля? Похоже, что вы совсем не ослабляете его в своем примере. Поскольку вы реализуете фильтр верхних частот, вам также необходимо установить нулевое значение постоянного тока.
Это объяснило бы низкочастотные искажения. У вас будет много пульсаций в частотной характеристике на низких частотах, если это значение постоянного тока будет отличным от нуля из-за большого перехода.
Вот пример в 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');
Обратите внимание, что в моем коде я создаю пример значения DC, отличного от нуля, за которым следует резкое изменение на ноль, и затем подъем. Затем я использую IFFT для преобразования во временную область. Затем я выполняю заполненный нулями fft (который автоматически выполняется MATLAB, когда вы передаете размер fft, больший, чем входной сигнал) для этого сигнала временной области. Заполнение нулями во временной области приводит к интерполяции в частотной области. Используя это, мы можем увидеть, как фильтр будет реагировать между выборками фильтров.
Одна из самых важных вещей, которую следует помнить, заключается в том, что даже если вы устанавливаете значения отклика фильтра на заданных частотах путем ослабления выходных сигналов ДПФ, это ничего не гарантирует для частот, возникающих между точками выборки.Это означает, что чем резче ваши изменения, тем больше будет перерегулирования и колебаний между сэмплами.
Теперь отвечу на ваш вопрос о том, как должна выполняться эта фильтрация. Есть несколько способов, но один из самых простых для реализации и понимания - это метод оконного дизайна. Проблема с вашим текущим дизайном в том, что ширина перехода огромна. В большинстве случаев вам потребуется как можно более быстрые переходы с минимально возможной рябью.
В следующем коде я создам идеальный фильтр и покажу ответ:
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');
Обратите внимание, что резкие изменения вызывают сильные колебания.
БПФ или дискретное преобразование Фурье - это дискретная версия преобразования Фурье. Преобразование Фурье применяется к сигналу в непрерывном диапазоне от бесконечности до бесконечности, в то время как ДПФ применяется к конечному числу выборок. Это, в сущности, приводит к квадратному окну (усечению) во временной области при использовании ДПФ, поскольку мы имеем дело только с конечным числом выборок. К сожалению, ДПФ прямоугольной волны является функцией типа sinc (sin (x) / x).
Проблема с резкими переходами в вашем фильтре (быстрый переход от 0 к 1 в одном сэмпле) заключается в том, что у него очень длинный отклик во временной области, которая усекается квадратным окном. Поэтому, чтобы минимизировать эту проблему, мы можем умножить сигнал во временной области на более постепенное окно.Если мы умножим окно обработки, добавив строку:
x = x .* hanning(1,N).';
после выполнения ОБПФ, мы получим такой ответ:
Поэтому я бы рекомендовал попробовать реализовать метод проектирования окон, поскольку он довольно прост (есть способы получше, но они усложняются). Поскольку вы реализуете эквалайзер, я предполагаю, что вы хотите иметь возможность изменять затухание на лету, поэтому я бы предложил рассчитывать и сохранять фильтр в частотной области всякий раз, когда происходит изменение параметров, а затем вы можете просто применить его. для каждого входного аудиобуфера, взяв fft входного буфера, умножив его на выборки фильтра частотной области, а затем выполнив ifft, чтобы вернуться во временную область. Это будет намного эффективнее, чем все ветвления, которые вы делаете для каждого образца.
Ваша основная проблема заключается в том, что частоты плохо определены в короткие промежутки времени . Это особенно верно для низких частот, поэтому вы больше всего замечаете проблему именно там.
Следовательно, когда вы берете действительно короткие сегменты из звуковой последовательности, а затем фильтруете их, отфильтрованные сегменты не будут фильтровать таким образом, чтобы получить непрерывную форму волны, и вы слышите скачки между сегментами, и это то, что генерирует щелкает вы здесь.
Например, возьмем некоторые разумные числа: я начинаю с формы волны с частотой 27,5 Гц (A0 на фортепиано), оцифрованной с частотой 44100 Гц, она будет выглядеть следующим образом (где красная часть имеет длину 1024 сэмпла):
Итак, сначала мы начнем с нижних частот 40 Гц. Итак, поскольку исходная частота меньше 40 Гц, фильтр нижних частот с отсечкой 40 Гц не должен иметь никакого эффекта, и мы получим выходной сигнал, который почти точно соответствует входному. Верно? Неправильно, неверно, неверно - и в этом суть вашей проблемы. Проблема в том, что для коротких участков идея 27,5 Гц четко не определена и не может быть хорошо представлена в ДПФ.
То, что 27,5 Гц не имеет особого значения в коротком сегменте, можно увидеть, взглянув на ДПФ на рисунке ниже. Обратите внимание, что хотя ДПФ более длинного сегмента (черные точки) показывает пик на частоте 27,5 Гц, короткий (красные точки) - нет.
Очевидно, что тогда фильтрация ниже 40 Гц будет просто фиксировать смещение постоянного тока, и результат фильтра нижних частот 40 Гц показан ниже зеленым цветом.
Синяя кривая (снятая с отсечкой 200 Гц) начинает соответствовать намного лучше.Но обратите внимание, что это не низкие частоты, которые заставляют его хорошо сочетаться, а включение высоких частот. Только когда мы включим все возможные частоты в короткий сегмент, вплоть до 22 кГц, мы наконец получим хорошее представление исходной синусоидальной волны.
Причина всего этого в том, что небольшой сегмент синусоидальной волны 27,5 Гц является не синусоидальной волной 27,5 Гц, и это DFT не имеет особого отношения к 27,5 Гц.