Как фильтровать данные FFT (для визуализации аудио)?

Если вы посмотрите на демонстрацию, пики fft падают плавно. Я пытаюсь сделать то же самое с обработкой в режиме Java, используя библиотеку minim. Я посмотрел, как это делается с помощью веб-аудио api в doFFTAnalysis() и попытался реплицировать это с помощью minim. Я также попытался передать, как abs() работает со сложным типом:

В настоящее время я делаю быстрый прототип с использованием Processing (java framework/library). Мой код выглядит следующим образом:

Я не получаю ошибок, что хорошо, но я не вижу ожидаемого поведения, которое плохо. Я ожидал, что пики будут падать медленнее, если сглаживание (k) близко к 1, но насколько я могу сказать только мой код масштабирует полосы.

К сожалению, математика и звук - это не моя сильная сторона, поэтому я ношу в темноте. Как я могу воспроизвести приятную визуализацию из демонстрации веб-аудио API?

У меня возникнет соблазн сказать, что это может быть языковым агностиком, но использование javascript, например, не применимо:). Тем не менее, я рад попробовать любую другую библиотеку java, которая выполняет анализ FFT.

У меня есть простое решение для сглаживания (непрерывно уменьшайте значения каждого предыдущего fft-диапазона, если текущий fft-диапазон не выше:

Увядший граф - сглаженный, а полностью насыщенный - живой.

Тем не менее, я все еще что-то пропущу, по сравнению с демонстрацией веб-аудио API:

Я думаю, что API веб-аудио может принимать во внимание, что средние и более высокие частоты необходимо масштабировать, чтобы быть ближе к тому, что мы воспринимаем, но я не уверен, как это решить.

Я пытался больше узнать о том, как класс RealtimeAnalyser делает это для WebAudioAPI, но кажется класс FFTFrame doFFT метод может сделать логарифмическое масштабирование. Я еще не понял, как работает makeFFT.

Как я могу масштабировать необработанный график FFT с логарифмической шкалой для учета восприятия? Моя цель - сделать приличную визуализацию, и я предполагаю, что мне нужно:

Любые подсказки о том, как я могу это достичь?

Я предполагаю, что эта часть выполняет сглаживание и масштабирование, которое я получаю в API веб-аудио: // Нормализуем так, как входная синусоидальная волна при 0dBfs регистрируется как 0dBfs (отменить коэффициент масштабирования FFT). const double величинаScale = 1.0/DefaultFFTSize;

Кажется, что масштабирование выполняется путем принятия абсолютного значения комплексного значения. Этот пост указывает в том же направлении. Я попытался использовать абс комплексного числа с помощью Minim и использовать различные функции окна, но он по-прежнему не похож на то, к чему я стремлюсь (Web Audio API):

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

Я попытался применить полезные советы @wakjah, например:

Я не уверен, что я применил подсказки по своему усмотрению. Вот как выглядит мой вывод:

но я не думаю, что я еще там, если я сравню это с визуализацией, к которой я стремился:

Я не уверен, что я правильно применил масштаб журнала. Мои предположения состояли в том, что я бы сделал сюжет, похожий на то, к чему я стремился после использования log10 (игнорируя сглаживание на данный момент).

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

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

Ответы

Ответ 1

Я немного неясен о том, какое именно сглаживание вы хотите сделать, но я попытаюсь предоставить некоторую информацию, которая может вам помочь.

Масштабирование результатов FFT для отображения

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

Чтобы сделать это масштабирование, мы преобразуем значения в децибелы. Важно отметить, что децибелы представляют собой шкалу, а не единицу - она представляет собой соотношение между двумя числами: обычно это измеренное значение и некоторая ссылка. Общая формула для децибел составляет

x_dB = 10 * log10((x ^ 2) / (ref ^ 2))

где log10 это логарифм по основанию 10, ^ является оператором мощности, а x_ref это выбранный опорное значение. Поскольку значения FFT'd из аудиофайла не имеют (обычно) каких-либо значимых единиц, x_ref обычно выбирается просто 1 для этого приложения. Далее, поскольку x является сложным, вам нужно принять абсолютное значение. Таким образом, формула будет

x_dB = 10 * log10(abs(x) ^ 2)

Здесь возможно небольшая (численная и скорость) оптимизация, поскольку вы возводите квадрат в квадратный корень:

x_dB = 10 * log10(real(x) ^ 2 + imag(x) ^ 2)

Перцептивное взвешивание

Масштабирование измерений частотной области обычно выполняется при измерении уровня звукового давления и мощности: для данного приложения выбран конкретный тип измерения (я не буду входить в типы здесь), и запись звука производится в соответствии с к этому типу измерения. Результатом является FFT'd, а затем умножается на заданное взвешивание на каждой частоте в зависимости от того, для чего будет использоваться результат и какой тип звука был записан. В общем случае используются два веса: A и C. C обычно используются только для звуков с чрезвычайно высокой амплитудой.

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

Информация о A-weighting в wikipedia.

Оконная

Окно выполняется в основном для уменьшения эффекта явления в этой бесплатной онлайн-книге.

Сглаживание во временной области отдельных ящиков частоты

Что касается того, что линия в каждом частотном блоке распадается медленно, вот простая идея, которая могла бы сделать трюк: в каждом частотном буфере применить простую экспоненциальную скользящую среднюю. Скажем, ваши результаты FFT хранятся в X[k], где k - это индекс частоты. Пусть ваше значение отображения Y[k] такое, что

Y[k] = alpha * Y_(t-1)[k] + (1 - alpha) * X[k]

где 0 < alpha < 1 - ваш коэффициент сглаживания, а Y_(t-1)[k] - значение Y[k] на последнем шаге времени (t-1). Это на самом деле простой фильтр низких частот IIR (бесконечный импульсный отклик), и, надеюсь, он должен делать в основном то, что вы хотите (возможно, с небольшой настройкой). Чем ближе альфа к нулю, тем быстрее будут влиять новые наблюдения (вход X[k]). Чем ближе он к одному, тем медленнее результат будет распадаться, но вход будет также влиять на результат медленнее, поэтому он может выглядеть "вялым". Вы можете захотеть добавить к нему условное условие, чтобы принять новое значение немедленно, если оно больше текущего значения.

Обратите внимание, что, опять же, это должно быть выполнено до преобразования в децибелы.

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

(править 2) Ваше третье обновление, опять же, близко, но есть небольшое неправильное преобразование формулы в следующих строках

fftSmooth[i] = smoothing * fftPrev[i] + ((1 - smoothing) * fftCurr[i]);

fftPrev[i] = fftCurr[i];//

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

fftSmooth[i] = smoothing * fftSmooth[i] + ((1 - smoothing) * fftCurr[i]);

Если smoothing == 0, эта строка должна иметь мало эффекта, кроме как умножить результат на скаляр.

Нормализация вычислений по абсолютной величине

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

Чтобы сделать это просто в коде, вы можете сделать что-то вроде

float maxVal = Math.max(Math.abs(fftReal[i]), Math.abs(fftImag[i]));
if (maxVal != 0.0f) { // prevent divide-by-zero
    // Normalize
    fftReal[i] = fftReal[i] / maxVal;
    fftImag[i] = fftImag[i] / maxVal;
}

fftCurr[i] = scale * (float)Math.log10(fftReal[i]*fftReal[i] + fftImag[i]*fftImag[i]);
// ...

Объединяя все это: Некоторый код

Некоторое время перепутавшись с ним в обработке 2.1, у меня есть решение, которое, я считаю, вам понравится:

import ddf.minim.analysis.*;
import ddf.minim.*;

Minim       minim;
//AudioInput  in;
AudioPlayer in;
FFT         fft;

float smoothing = 0.60;
final boolean useDB = true;
final int minBandwidthPerOctave = 200;
final int bandsPerOctave = 10;
float[] fftSmooth;
int avgSize;

float minVal = 0.0;
float maxVal = 0.0;
boolean firstMinDone = false;

void setup(){
  minim = new Minim(this);
  //in = minim.getLineIn(Minim.STEREO, 512);
  in = minim.loadFile("C:\\path\\to\\some\\audio\\file.ext", 2048);

  in.loop();

  fft = new FFT(in.bufferSize(), in.sampleRate());

  // Use logarithmically-spaced averaging
  fft.logAverages(minBandwidthPerOctave, bandsPerOctave);

  avgSize = fft.avgSize();
  fftSmooth = new float[avgSize];

  int myWidth = 500;
  int myHeight = 250;
  size(myWidth, myHeight);
  colorMode(HSB,avgSize,100,100);

}

float dB(float x) {
  if (x == 0) {
    return 0;
  }
  else {
    return 10 * (float)Math.log10(x);
  }
}

void draw(){
  background(0);
  stroke(255);

  fft.forward( in.mix);

  final int weight = width / avgSize;
  final float maxHeight = (height / 2) * 0.75;

  for (int i = 0; i < avgSize; i++) {
    // Get spectrum value (using dB conversion or not, as desired)
    float fftCurr;
    if (useDB) {
      fftCurr = dB(fft.getAvg(i));
    }
    else {
      fftCurr = fft.getAvg(i);
    }

    // Smooth using exponential moving average
    fftSmooth[i] = (smoothing) * fftSmooth[i] + ((1 - smoothing) * fftCurr);

    // Find max and min values ever displayed across whole spectrum
    if (fftSmooth[i] > maxVal) {
      maxVal = fftSmooth[i];
    }
    if (!firstMinDone || (fftSmooth[i] < minVal)) {
      minVal = fftSmooth[i];
    }
  }

  // Calculate the total range of smoothed spectrum; this will be used to scale all values to range 0...1
  final float range = maxVal - minVal;
  final float scaleFactor = range + 0.00001; // avoid div. by zero

  for(int i = 0; i < avgSize; i++)
  {
    stroke(i,100,100);
    strokeWeight(weight);

    // Y-coord of display line; fftSmooth is scaled to range 0...1; this is then multiplied by maxHeight
    // to make it within display port range
    float fftSmoothDisplay = maxHeight * ((fftSmooth[i] - minVal) / scaleFactor);

    // X-coord of display line
    float x = i * weight;

    line(x, height / 2, x, height / 2 - fftSmoothDisplay);
  }
  text("smoothing: " + (int)(smoothing*100)+"\n",10,10);
}
void keyPressed(){
  float inc = 0.01;
  if(keyCode == UP && smoothing < 1-inc) smoothing += inc;
  if(keyCode == DOWN && smoothing > inc) smoothing -= inc;
}

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

Улучшение: теперь с A-weighting

У меня есть обновленная версия кода, который применяет взвешивание A в каждом частотном диапазоне (хотя только в том случае, когда режим dB включен, потому что таблица была у меня в дБ:). Поверните A-weighting для получения результата ближе к WMP или выключите его ближе к VLC.

Есть также некоторые незначительные изменения в том, как он отображается: он теперь центрирован на дисплее, и он отображает только максимальную центральную частоту диапазона.

Здесь код - наслаждайтесь!

import ddf.minim.analysis.*;
import ddf.minim.*;

Minim       minim;
//AudioInput  in;
AudioPlayer in;
FFT         fft;

float smoothing = 0.73;
final boolean useDB = true;
final boolean useAWeighting = true; // only used in dB mode, because the table I found was in dB 
final boolean resetBoundsAtEachStep = false;
final float maxViewportUsage = 0.85;
final int minBandwidthPerOctave = 200;
final int bandsPerOctave = 10;
final float maxCentreFrequency = 18000;
float[] fftSmooth;
int avgSize;

float minVal = 0.0;
float maxVal = 0.0;
boolean firstMinDone = false;

final float[] aWeightFrequency = { 
  10, 12.5, 16, 20, 
  25, 31.5, 40, 50, 
  63, 80, 100, 125, 
  160, 200, 250, 315, 
  400, 500, 630, 800, 
  1000, 1250, 1600, 2000, 
  2500, 3150, 4000, 5000,
  6300, 8000, 10000, 12500, 
  16000, 20000 
};

final float[] aWeightDecibels = {
  -70.4, -63.4, -56.7, -50.5, 
  -44.7, -39.4, -34.6, -30.2, 
  -26.2, -22.5, -19.1, -16.1, 
  -13.4, -10.9, -8.6, -6.6, 
  -4.8, -3.2, -1.9, -0.8, 
  0.0, 0.6, 1.0, 1.2, 
  1.3, 1.2, 1.0, 0.5, 
  -0.1, -1.1, -2.5, -4.3, 
  -6.6, -9.3 
};

float[] aWeightDBAtBandCentreFreqs;

void setup(){
  minim = new Minim(this);
  //in = minim.getLineIn(Minim.STEREO, 512);
  in = minim.loadFile("D:\\Music\\Arthur Brown\\The Crazy World Of Arthur Brown\\1-09 Fire.mp3", 2048);

  in.loop();

  fft = new FFT(in.bufferSize(), in.sampleRate());

  // Use logarithmically-spaced averaging
  fft.logAverages(minBandwidthPerOctave, bandsPerOctave);
  aWeightDBAtBandCentreFreqs = calculateAWeightingDBForFFTAverages(fft);

  avgSize = fft.avgSize();
  // Only use freqs up to maxCentreFrequency - ones above this may have
  // values too small that will skew our range calculation for all time
  while (fft.getAverageCenterFrequency(avgSize-1) > maxCentreFrequency) {
    avgSize--;
  }

  fftSmooth = new float[avgSize];

  int myWidth = 500;
  int myHeight = 250;
  size(myWidth, myHeight);
  colorMode(HSB,avgSize,100,100);

}

float[] calculateAWeightingDBForFFTAverages(FFT fft) {
  float[] result = new float[fft.avgSize()];
  for (int i = 0; i < result.length; i++) {
    result[i] = calculateAWeightingDBAtFrequency(fft.getAverageCenterFrequency(i));
  }
  return result;    
}

float calculateAWeightingDBAtFrequency(float frequency) {
  return linterp(aWeightFrequency, aWeightDecibels, frequency);    
}

float dB(float x) {
  if (x == 0) {
    return 0;
  }
  else {
    return 10 * (float)Math.log10(x);
  }
}

float linterp(float[] x, float[] y, float xx) {
  assert(x.length > 1);
  assert(x.length == y.length);

  float result = 0.0;
  boolean found = false;

  if (x[0] > xx) {
    result = y[0];
    found = true;
  }

  if (!found) {
    for (int i = 1; i < x.length; i++) {
      if (x[i] > xx) {
        result = y[i-1] + ((xx - x[i-1]) / (x[i] - x[i-1])) * (y[i] - y[i-1]);
        found = true;
        break;
      }
    }
  }

  if (!found) {
    result = y[y.length-1];
  }

  return result;     
}

void draw(){
  background(0);
  stroke(255);

  fft.forward( in.mix);

  final int weight = width / avgSize;
  final float maxHeight = height * maxViewportUsage;
  final float xOffset = weight / 2 + (width - avgSize * weight) / 2;

  if (resetBoundsAtEachStep) {
    minVal = 0.0;
    maxVal = 0.0;
    firstMinDone = false;
  }

  for (int i = 0; i < avgSize; i++) {
    // Get spectrum value (using dB conversion or not, as desired)
    float fftCurr;
    if (useDB) {
      fftCurr = dB(fft.getAvg(i));
      if (useAWeighting) {
        fftCurr += aWeightDBAtBandCentreFreqs[i];
      }
    }
    else {
      fftCurr = fft.getAvg(i);
    }

    // Smooth using exponential moving average
    fftSmooth[i] = (smoothing) * fftSmooth[i] + ((1 - smoothing) * fftCurr);

    // Find max and min values ever displayed across whole spectrum
    if (fftSmooth[i] > maxVal) {
      maxVal = fftSmooth[i];
    }
    if (!firstMinDone || (fftSmooth[i] < minVal)) {
      minVal = fftSmooth[i];
    }
  }

  // Calculate the total range of smoothed spectrum; this will be used to scale all values to range 0...1
  final float range = maxVal - minVal;
  final float scaleFactor = range + 0.00001; // avoid div. by zero

  for(int i = 0; i < avgSize; i++)
  {
    stroke(i,100,100);
    strokeWeight(weight);

    // Y-coord of display line; fftSmooth is scaled to range 0...1; this is then multiplied by maxHeight
    // to make it within display port range
    float fftSmoothDisplay = maxHeight * ((fftSmooth[i] - minVal) / scaleFactor);
    // Artificially impose a minimum of zero (this is mathematically bogus, but whatever)
    fftSmoothDisplay = max(0.0, fftSmoothDisplay);

    // X-coord of display line
    float x = xOffset + i * weight;

    line(x, height, x, height - fftSmoothDisplay);
  }
  text("smoothing: " + (int)(smoothing*100)+"\n",10,10);
}
void keyPressed(){
  float inc = 0.01;
  if(keyCode == UP && smoothing < 1-inc) smoothing += inc;
  if(keyCode == DOWN && smoothing > inc) smoothing -= inc;
}

Ответ 2

Внутри вашего цикла: вам нужно добавить логарифмическое вычисление для шкалы lg:

stroke(i,100,50);
line( i, height, i, height - fftSmooth[i]*8 );
stroke(i,100,100);
line( i, height, i, height - band*8 );

Должно быть изменено на:

int l = map(log(map(i ,0 ,specSize,0,100),0,2,0,width).  // an estimation, may have to calibrate
stroke(i,100,50);
line( l, height, l, height - fftSmooth[i]*8 );
stroke(i,100,100);
line( l, height, l, height - band*8 );