Надо получить частотну характеристику звука. Как я делаю ниже.
Проблема в том, что в качестве пиков выдаются какие-то ну очень приблизительные значения (разброс 100-200 герц), хотя, скажем, soundforge по тому же сигналу определяет пики очень чётко.
Формат записи 8бит, 22500Гц, моно.
static void dsAudioCapture_OnCapture(CaptureParameters captureParameters, byte[] data)
{
/* получили новые данные и добавили в буффер */foreach (byte sample in data)
{
buffer.Add(sample);
}
/* если данных в буффере достаточно для обработки */if (buffer.Count > BlockSize)
{
/* найдём минимальное и максимальное значения */double bufferMin = buffer[0];
double bufferMax = buffer[0];
for (int index = 0; index < BlockSize; ++index)
{
if (buffer[index] < bufferMin)
{
bufferMin = buffer[index];
}
if (buffer[index] > bufferMax)
{
bufferMax = buffer[index];
}
}
/* тут хранится сам сигнал */
Complex[] signal = new Complex[BlockSize];
/* отмасштабируем данные, чтобы работать с числами определённого порядка и избежать ошибок округления */for (int index = 0; index < BlockSize; ++index)
{
signal[index] = new Complex((buffer[index] - bufferMin) / (bufferMax - bufferMin), 0.0);
}
/* выкинем обработанный блок из буффера */
buffer.RemoveRange(0, BlockSize);
/* тут будут записаны результаты */
Complex[] result = new Complex[BlockSize];
/* натравили БПФ */
fft.Calculate(signal, result);
/* тут будут храниться пары интенсивность-частота */
SortedList<double, double> frequencyList = new SortedList<double, double>();
/* интересуют частоты ниже 2КГц */int indexTo = (int)Math.Round(2000.0 * BlockSize / captureParameters.SamplePerSecondNumber);
/* заполняем сортированный список частот. Самые интенсивные в конце */for (int index = 1; index < indexTo; ++index)
{
frequencyList.Add(
result[index].Re,
((double)index) * ((double)captureParameters.SamplePerSecondNumber) / ((double)BlockSize));
}
/* только интересующие частоты */int frequencyFrom = 500;
int frequencyTo = 2000;
int frequencyCount = 0;
/* есть ли среди 5 самых интенсивных частот интересующие? */for (int index = Math.Max(0, frequencyList.Count - 5); index < frequencyList.Count; ++index)
{
if ((frequencyList.Values[index] > frequencyFrom) && (frequencyList.Values[index] < frequencyTo))
{
++frequencyCount;
}
}
/* если да, то выводим их */if (frequencyCount > 0)
{
Console.Write("FREQ:");
for (int index = Math.Max(0, frequencyList.Count - 5); index < frequencyList.Count; ++index)
{
if ((frequencyList.Values[index] > frequencyFrom) && (frequencyList.Values[index] < frequencyTo))
{
Console.Write(" [{0}] ", Math.Round(frequencyList.Values[index]).ToString().PadLeft(5));
}
}
Console.WriteLine();
}
}
}
Во избежание лишних вопросов, привожду так же код класса БПФ (Совершенно не оптимизирован, но мне пока важно, чтобы всё работало. E — even, O — odd.).
class FFT
{
public void Calculate(Complex[] signal, Complex[] result)
{
if (signal.Length == 1)
{
result[0] = signal[0];
}
else
{
Complex[] signalE = new Complex[signal.Length / 2];
Complex[] signalO = new Complex[signal.Length / 2];
for (int index = 0; index < signal.Length / 2; ++index)
{
signalE[index] = signal[2 * index];
signalO[index] = signal[2 * index + 1];
}
Complex[] resultE = new Complex[signal.Length / 2];
Complex[] resultO = new Complex[signal.Length / 2];
Calculate(signalE, resultE);
Calculate(signalO, resultO);
Complex omegaN = new Complex(Math.Cos(2 * Math.PI / signal.Length), Math.Sin(2 * Math.PI / signal.Length));
Complex omega = new Complex(1.0);
for (int index = 0; index < signal.Length / 2; ++index)
{
result[index] = resultE[index] + resultO[index] * omega;
result[index + signal.Length / 2] = resultE[index] - resultO[index] * omega;
omega = omega * omegaN;
}
}
}
}
Ну то что я Complex реализовал правильно, надеюсь вы мне верите
adontz wrote: > Надо получить частотну характеристику звука. Как я делаю ниже. > Проблема в том, что в качестве пиков выдаются какие-то ну очень > приблизительные значения (разброс 100-200 герц), хотя, скажем, > soundforge по тому же сигналу определяет пики очень чётко. > Формат записи 8бит, 22500Гц, моно.
Стандартные грабли — вроде, окно не используется. Дело в том, что БПФ от
отрезка звука — это БПФ его свёртки с прямоугольником,а при этом спектр
умножается на (sin x)/x . Поэтому для заглаживания краёв перед БПФ
сигнал умножают на окно — функцию, плавно спадающую до нуля в концах
отрезка (желательно — значение и производная по нулям, кажется) и
имеющую быстро спадающий спектр.
В справке MatLab, http://www.mathworks.com/access/helpdesk/help/toolbox/signal/reftabl2.html#136851
приведены формулы для многих оконных функций.
В учебнике http://www.bores.com/courses/advanced/windows/ объяснено про
окна чуть подробнее.
Здравствуйте, adontz, Вы писали:
A>Богатый выбор иногда рождает больше проблем, чем отсутствие выбора A>Какую оконную функцию брать-то?
Что-то мне кажется, что сама задача как-то неверно поставлена. Что требуется-то? Определить частотные харктеристики на некоей регулярной сетке? Или выдавать мнгновенные (и нерегулярные, т.е. дробные) величины и фазы текущих частотных составляющих на каждом отсчете (скажем до 10 основных по можности тонов?)
Решал когда-то подобную задачу для электрогитарного синтезатора. Есть алгоритмы. Они обычно бесконечны, работают без окна и всяких БПФ, т.е. организованы как бесконечные фильтры, "следящие" за основными (по мощности) составляющими звука. Как раз мне под гитарный процессор подходило, ибо нужна была мгновенная точность основных составляющих.
Здравствуйте, adontz, Вы писали:
A>Здравствуйте, vdimas, Вы писали:
V>>Что-то мне кажется, что сама задача как-то неверно поставлена. Что требуется-то?
A>Распознавать DTMF
Это не делают БПФ ом. Пороюсь в своих архивах, напишу позже.
Здравствуйте, adontz, Вы писали:
A>Здравствуйте, vdimas, Вы писали:
V>>Что-то мне кажется, что сама задача как-то неверно поставлена. Что требуется-то?
A>Распознавать DTMF
Я так понимаю, что нужно определять основную частоту сигнала. Для этого обычно используются методы параметрического спектрального анализа (например, MUSIC или EV). Точность получаемых результатов у них больше. Про эти методы можно, например почитать в "Цифровой спектральный анализ и его приложения" (Марпл-мл.).
Здравствуйте, adontz, Вы писали:
A>Здравствуйте, jhng, Вы писали:
A>>>Распознавать DTMF J>>Я так понимаю, что нужно определять основную частоту сигнала.
A>Dual-Tone Multi-Frequency
A>Нужно найти интенсивности всех гармоник. Просто найти частоту основной гармоники мало.
Ну так ты прочитай про эти методы поподробнее. Они как раз и предназначены для определения частот и уровней гармонических составляющих сигнала.
Если все возможные частоты сигналов в канале зараннее известны, то можно использовать набор оптимальных фильтров (фильтров Винера) и анализировать уровень сигнала на выходе каждого из них. Затем по получаемым результатам принимать то или иное решение. Может так будет даже проще.
Здравствуйте, jhng, Вы писали:
J>Здравствуйте, adontz, Вы писали:
A>>Здравствуйте, jhng, Вы писали:
A>>>>Распознавать DTMF J>>>Я так понимаю, что нужно определять основную частоту сигнала.
A>>Dual-Tone Multi-Frequency
A>>Нужно найти интенсивности всех гармоник. Просто найти частоту основной гармоники мало.
J>Ну так ты прочитай про эти методы поподробнее. Они как раз и предназначены для определения частот и уровней гармонических составляющих сигнала. J>Если все возможные частоты сигналов в канале зараннее известны, то можно использовать набор оптимальных фильтров (фильтров Винера) и анализировать уровень сигнала на выходе каждого из них. Затем по получаемым результатам принимать то или иное решение. Может так будет даже проще.
А можно еще про вейвлеты почитать здесь, они как раз и предназначены для получения частотной информации по времени.
Здравствуйте, jhng, Вы писали:
J>Здравствуйте, adontz, Вы писали:
A>>Здравствуйте, jhng, Вы писали:
A>>>>Распознавать DTMF J>>>Я так понимаю, что нужно определять основную частоту сигнала.
A>>Dual-Tone Multi-Frequency
A>>Нужно найти интенсивности всех гармоник. Просто найти частоту основной гармоники мало.
Еще как вариант можно использовать банк резонансных контуров (с резистором в цепи конденсатора или катушки для ограничения коэффициента передачи по частоте). Тогда и групповая задержка у всех фильтров будет одинаковой на резонансной частоте (наверное... ), и затрат ресурсов/времени мало, поскольку разностные уравнения получаются простыми.
Здравствуйте, adontz, Вы писали:
A>Здравствуйте, Шахтер, Вы писали:
Ш>>Это не делают БПФ ом. Пороюсь в своих архивах, напишу позже.
A>Да? А у меня почти вышло
Держи код (не тестировал).
Если коротко -- то мы считаем энергию синусоидальной составляющей данной частоты (по 8 частотам). При поступлении нового отсчета надо сделать небольшой пересчет. Алгоритм получается очень экономным.
namespace DTMFDetect {
struct DTMFCode
{
int low; // 1-4 , 0 int hi; // 1-4 , 0 public DTMFCode(bool fake/*false*/) { low=hi=0; }
public DTMFCode(int low_,int hi_) { low=low_; hi=hi_; }
}
class DTMFDetector
{
const int FreqCount=8;
const int BufSize=512; // Число отчетов фрейма сигнала, должно быть достаточным, чтобы вмещать
// 10 или более периодов синуса искомых частот.
// Увеличение этого числа приводит к увеличению селективности детектора. const double Fd=22000; // Частота дискретизации.const double Scale=(1<<15)-2;
const double ELim=BufSize; // отсечение по энергии const double ProdLim=((Scale*Scale*BufSize)/4)*0.9; // чувствительность по частоте
Int16[] buf;
int index;
bool full;
Int16[,] sinTable;
Int16[,] cosTable;
Int64 E;
Int64[] sinProd;
Int64[] cosProd;
static double Norm(double x,double y) { return x*x+y*y; }
DTMFCode check()
{
if( E>ELim )
{
for(int i=0; i<4 ;i++)
{
if( Norm(sinProd[i],cosProd[i])>E*ProdLim )
{
for(int j=4; j<8 ;j++)
{
if( Norm(sinProd[j],cosProd[j])>E*ProdLim )
{
return new DTMFCode(i+1,j-4+1);
}
}
}
}
}
return new DTMFCode(false);
}
public DTMFDetector()
{
buf=new Int16[BufSize];
sinTable=new Int16[FreqCount,BufSize];
cosTable=new Int16[FreqCount,BufSize];
index=0;
full=false;
E=0;
sinProd=new Int64[FreqCount];
cosProd=new Int64[FreqCount];
int i=0;
foreach(double F in new double[]{ 697,770,852,941, 1209,1336,1477,1633 })
{
double omega=2*Math.PI*(F/Fd);
for(int j=0; j<BufSize ;j++)
{
sinTable[i,j]=(Int16)( Scale*Math.Sin(omega*j) );
cosTable[i,j]=(Int16)( Scale*Math.Cos(omega*j) );
}
i++;
}
}
public DTMFCode next(Int16 x)
{
Int16 y=buf[index];
//
// Фрейм сигнала хранится в циклическом буфере
// При поступлении отсчета пересчитываем энергию фрейма и скалярные произведения с синусами и косинусами
//
E+=((Int64)x*x-(Int64)y*y);
for(int i=0; i<FreqCount ;i++)
{
sinProd[i]+=((Int64)x-y)*sinTable[i,index];
cosProd[i]+=((Int64)x-y)*cosTable[i,index];
}
buf[index]=x;
if( ++index>=BufSize )
{
index=0;
full=true;
}
return full?check():new DTMFCode(false);
}
}
} // namespace DTMFDetect
J>Если все возможные частоты сигналов в канале зараннее известны, то можно использовать набор оптимальных фильтров (фильтров Винера) и анализировать уровень сигнала на выходе каждого из них.
оптимальный фильтр для синусоиды как выглядит?
Здравствуйте, LelicDsp, Вы писали:
J>>Если все возможные частоты сигналов в канале зараннее известны, то можно использовать набор оптимальных фильтров (фильтров Винера) и анализировать уровень сигнала на выходе каждого из них. LD>оптимальный фильтр для синусоиды как выглядит?
Берешь один период синусоиды, вычисляешь ее АКФ. Полученные отсчеты АКФ и будут являться коэффициентами фильтра.
Нечто похожее Шахтер уже привел.
Здравствуйте, adontz, Вы писали:
A>Здравствуйте, jhng, Вы писали:
J>>Берешь один период синусоиды, вычисляешь ее АКФ
A>Что такое АКФ? Люди я в этом вопросе лопух, не надо кидаться непонятными словами
Во первых sorry за неправильную информацию в предыдущем посте. Оптимальным для выбранного сигнала будет являться такой фильтр, коэффициенты которого будут являться зеркальным отражением отсчетов исходного сигнала за период. Так что для грамонического сигнала sin(omega * t + phi) коэффициенты фильтра можно расчитать как 1 — sin(omega * t + phi). Соответственно сигнал на выходе будет совпадать с АКФ входного сигнала. В приведенный выше коде и реализован этот подход. Но все таки мне использование резонаторов кажеться лучше, поскольку они проще в реализации и могут дать лучшее разрешение по чатоте (а частоты в DTMF расположены очень близко). Однако это только в том случае, если отношение сигнал/шум на частоте резонанса достаточно велико (скорее всего это условие выполняется).
Собственно про АКФ. АКФ — автокорреляционная функция. Она показывает степень сходства сигнала со своей сдвинутой копией. Расчитывается по формуле R(t) = int(s(t) * s(t — tau) dtau). Пределы интегрирования [-inf; +inf]. В дискретной области АКФ расчитывается аналогично, только интеграл заменяется на сумму.
P.S. Поковырял в matlab методы спектрального оценивания. Без шума частоты оценивает просто супер.
А вот если шум добавить, тогда... Короче если даже небольшой шум в канале связи (-20 dB) — не прокатит, слишком близко частоты.
P.P.S. Если с математикой не очень — скажи все частоты DTMF и частоту дискретизации. Может сам фильтры посчитаю, если время будет
Здравствуйте, adontz, Вы писали:
V>>Что-то мне кажется, что сама задача как-то неверно поставлена. Что требуется-то?
A>Распознавать DTMF
Это делается коррелляцией, даже точность синусоид никакая не нужна, распознает железно на корреляции с прямоугольными импульсами. У меня DTMF распознавали старенькие 8039 и одновременно делали кучу других вещей.
В общем, на каждую частоту из DTMF заводится 2 фильтра НЧ первого порядка (т.е. на каждом шаге что-то типа s=s-(s/4) ) и 2 ячейки-сумматора.
Делаем непрерывную корреляцию с прямоугольными импульсами, сдвинутыми на 90 градусов для каждой частоты. Т.е. просто прибавляем к сумме (если "1") или отнимаем ("0"). Подсчет накопительных сумм корреляций обычно производится для каждой частоты индивидуально, кратно целому числу периодов (1 период — нормально, но не обязательно). Затем эти суммы фильтруются упомянутыми фильтрами.
Искомая текущая мощность сигнала на каждой из исследуемых частот — это корень квадратный по теореме пифагора: sqrt(abs(ФНЧ1)^2+abs(ФНЧ2)^2). Считать корень квадратный вовсе не обязательно. Достаточно работать непосредственно с квадратами мощностей. Задача сводится к подбору порога отношения мощностей исследуемых частот к общей мощности входного сигнала. Вернее, к подбору квадрата этого отношения. (Разумеется, можно так же ограничить абсолютным порогом "чуствительности" снизу всю систему)
Как генерировать прямоугольные импульсы, сдвинутые на 90 градусов.
предположим у тебя 8кГц тактовая, а тебе нужна частота 1200Гц. Берем накопитель 2 байта, например, т.е. 65536. Считаем: 1200*65536/8000 = 9830,4 (грубо — 9830, для нашей задачи точности хватает, если у тебя 4-байтное слово, то точность получается просто астрономической).
Т.е. прибавляешь 9830 к накопителю каждый раз, старший бит — это искомая частота 1200. Чтобы получить сдвиг на 90 градусов — это XOR со следующим битом.
