[rudy_b 1]

Интересно, спасибо. А еще, мне всегда казалось, что струна выдает не чистую синусоиду, а все-таки с гармониками, как минимум с четными, то есть у струны отдельно есть колебания ее половины, четверти...

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

 

Т.е. если есть моды с частотами 1 кГц и 1.001 кГц, то погрешность определения частоты может достигать 1 Гц. Даже формально в этом случае непонятно, что называть резонансной частотой датчика. А подобное наблюдалось - неправильное крепление струны приводит к появлению сильно разнесенных по частоте мод. Спасает только селективное возбуждение только одной из них - а это уже не ударное возбуждение, а долгая (более 1 сек) накачка соответствующей частотой.

[tmtlib 0]

Если кому ещё интересна эта тема, вот хороший материальчик:

http://stackoverflow.com/questions/4633203...-between-frames

 

Метод уточнения частоты с использованием данных о фазе (нужно брать 2 фурье, лучше с нахлёстом):

const double freqPerBin = SAMPLE_RATE / FFT_N;
const double phaseStep = 2.0 * M_PI * FFT_STEP / FFT_N;

// process phase difference
double delta = phase - m_fftLastPhase[k];
m_fftLastPhase[k] = phase;
delta -= k * phaseStep;  // subtract expected phase difference
delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval
delta /= phaseStep;  // calculate diff from bin center frequency
double freq = (k + delta) * freqPerBin;  // calculate the true frequency

 

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

[fontp 0]

Если кому ещё интересна эта тема, вот хороший материальчик:

http://stackoverflow.com/questions/4633203...-between-frames

 

Метод уточнения частоты с использованием данных о фазе (нужно брать 2 фурье, лучше с нахлёстом):

const double freqPerBin = SAMPLE_RATE / FFT_N;
const double phaseStep = 2.0 * M_PI * FFT_STEP / FFT_N;

// process phase difference
double delta = phase - m_fftLastPhase[k];
m_fftLastPhase[k] = phase;
delta -= k * phaseStep;  // subtract expected phase difference
delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval
delta /= phaseStep;  // calculate diff from bin center frequency
double freq = (k + delta) * freqPerBin;  // calculate the true frequency

 

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

 

Хороший алгоритм по типу штангенциркуля. Только Вы забыли сказать, что такое k. Из-за этого, если не читать матерьяльчик, то непонятно, что здесь написано. Здесь k - это номер бина в котором максимум энергии. В матерьяльчике автор показывает, что правильную частоту можно получить не только по этому "максимальному" бину но и по FFT_N/FFT_STEP соседних. Поэтому не возникает проблем с частотами на границе между бинами при перекрытии блоков FFT. Без нахлёста с такими пограничными частотами алгоритм работать не будет, шумы не дадут. Систематическая ошибка у алгоритма отсутствует. Интересно посмотреть как алгоритм ведёт себя по отношению к шуму - судя по всему выйдет на точность оценки максимального правдоподобия CRLB, особенно если ещё и научиться делать взвешенную оценку по нескольким значениям возле максимума (остальные FFT_N/FFT_STEP оценок из-за малой энергии бинов реально утонут под шумом)

[ivan219 0]

Если кому ещё интересна эта тема, вот хороший материальчик:

http://stackoverflow.com/questions/4633203...-between-frames

 

Метод уточнения частоты с использованием данных о фазе (нужно брать 2 фурье, лучше с нахлёстом):

const double freqPerBin = SAMPLE_RATE / FFT_N;
const double phaseStep = 2.0 * M_PI * FFT_STEP / FFT_N;

// process phase difference
double delta = phase - m_fftLastPhase[k];
m_fftLastPhase[k] = phase;
delta -= k * phaseStep;  // subtract expected phase difference
delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval
delta /= phaseStep;  // calculate diff from bin center frequency
double freq = (k + delta) * freqPerBin;  // calculate the true frequency

 

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

А можно этот код расписать иначе? Так как в матлабе не силён. Если можно то блок схемой или просто словами.

[tmtlib 0]

А можно этот код расписать иначе? Так как в матлабе не силён. Если можно то блок схемой или просто словами.

Я планирую выложить примерчик с интересными методами здесь. Но думаю это будет не скоро. А пока попробую на словах, например, для фурье на 1024 отсчёта:

1) Берём FFT и получаем массивы амплитуд мнимой и действительной части по 512 элементов в каждом Im[0..511], Re[0..511].

2) Для каждого из 512 элемента считаем фазу. Если по простому, то фаза:=ArcTan(Im / Re); Можно сделать получше, расписав разные случаи (если Re=0 и т.д.)

3) храним фазы за предыдущий проход FFT и за текущий. m_fftLastPhase - это старая фаза, а phase - это на самом деле phase, т.е. фаза от текущего FFT. Заметил, что метод лучше работает если применить сглаживающее окно.

 

Забыл написать: например, константа SAMPLE_RATE=44100, FFT_N =1024, FFT_STEP = 100 (количество точек, на которые смещаем исходные данные. 100 намного меньше 1024, поэтому результат лучше)

Изменено 13 ноября, 2011 пользователем tmtlib

[ivan219 0]

И всё таки не много не понятно.

Что делает эта функция?

delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval

[thermit 1]

вычисляет остаток от деления дельты на два пи, вероятно...

Изменено 13 ноября, 2011 пользователем thermit

[tmtlib 0]

И всё таки не много не понятно.

Что делает эта функция?

delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval

 

Фазу приводят к диапазону -3.14...+3.14, пример в градусах: это что-то наподобие того, что 350 градусов привратится в -10 и т.д. Т.е. вектор (Im,Re) всегда в правом полукруге от -90 до 90 градусов.

 

[fontp 0]

Фазу приводят к диапазону -3.14...+3.14, пример в градусах: это что-то наподобие того, что 350 градусов привратится в -10 и т.д. Т.е. вектор (Im,Re) всегда в правом полукруге от -90 до 90 градусов.

 

Именно на этом этапе алгоритм перестанет работать, если нет перекрытия. Как показывает автор "матерьяльчика" в последней табличке (Pass 4) для отсутствия перекрытия.

http://www.dspdimension.com/admin/pitch-sh...g-using-the-ft/

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

 

2) Для каждого из 512 элемента считаем фазу. Если по простому, то фаза:=ArcTan(Im / Re); Можно сделать получше, расписав разные случаи (если Re=0 и т.д.)

3) храним фазы за предыдущий проход FFT и за текущий. m_fftLastPhase - это старая фаза, а phase - это на самом деле phase, т.е. фаза от текущего FFT. Заметил, что метод лучше работает если применить сглаживающее окно.

 

Для определения частоты вам не нужны все FFT_N значений фазы, поскольку только FFT_N/FFT_STEP значений фазы около максимума спектральной мощности правильные, а остальные искажены врэпингом на 2*pi (Pass 5). Более того в реальной ситуации присутствуют шумы и большинство из этих FFT_N/FFT_STEP правильных оценок (если их много) будут реально сильно зашумлены из-за очень быстрого падения уровня сигнала в бине при отклонении от центрального. Мы помним, что отклик бинов на синусоиду падает в соответствии с функцией спектрального окна, и в лучшем случае он будет спадать как синк или быстрее. В реальной ситуации только 2-4 центральных бина будут давать значимую оценку частоты, причем наименее искажены шумом будут только одна (если частота близко к центру бина) или две (если частота на границе бинов) оценки для центральных бинов (где максимум спектральной мощности). Остальные "утонут" под шумом

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

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

http://www.dspdimension.com/admin/pitch-sh...g-using-the-ft/

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

 

С учетом сказанного простейший алгоритм будет выглядеть так:

1) Берём FFT и получаем массивы амплитуд мнимой и действительной части по 512 элементов в каждом Im[0..511], Re[0..511].

2) Вычисляем энергию для каждого k E[k]=Re*Re+Im*Im

3) Находим максимум E[k] - пусть будет E() для kmax. Здесь вообще-то нужно проверять что kmax не прыгает. Если изменение kmax от фрейма к фрейму больше 1 по модулю - то это не синусоида, выход по ошибке

4) Для элементов kmax, kmax-1, kmax+1 считаем фазу. фаза[]:=ATAN2(Im[], Re[]); другие элементы на самом деле не нужны, в них отсутствует полезная информация

3) храним посчитаные фазы за предыдущий проход FFT и за текущий. m_fftLastPhase - это старая фаза, а phase - это на самом деле phase[kmax], т.е. фаза[kmax] от текущего FFT. Поскольку kmax гарантировано не меняется больше чем на 1, то в любом случае значение фазы для kmax посчитаны и для предыдущего фрейма и есть в массиве

 

Если перекрытие FFT_N/FFT_STEP >2 то оценку можно немного улучшить если использовать не одно значение максимума, а два первых максимума и построить взвешенную оценку для частоты по этим двум оценкам.

Если главный максимальный бин k0=kmax, а следующий по энергии k1=arg max(E(kmax-1), E(kmax+1)), то можно построить оценку типа

(f(k1)*E(k1)+f(k0)*E(k0))/(E(k0)+E(k1)) с большим иммунитетом к шуму, главным образом для частот на границе бинов (поскольку в этом случае будет две равноценных оценки)

 

А так да, идея алгоритма на вскидку смотрится очень хорошо

[alex_os 0]

...

(f(k1)*E(k1)+f(k0)*E(k0))/(E(k0)+E(k1)) с большим иммунитетом к шуму, главным образом для частот на границе бинов (поскольку в этом случае будет две равноценных оценки)

 

f(k1) это измеренная частота для бина k1 ? А что если приращение фазы измерять по нескольким бинам в окрестности

максимума? Взвешивание получится автоматически.

 

        delta =  angle( X_1(kmax-n:kmax+n)' * X(kmax-n:kmax+n) )

 

X_1 - предидущий выход FFT (вектор столбец),

X - текущий выход FFT,

kmax - индекс бина с максимальной энергией,

n - число 0...3

[fontp 0]

f(k1) это измеренная частота для бина k1 ? А что если приращение фазы измерять по нескольким бинам в окрестности

максимума? Взвешивание получится автоматически.

 

        delta =  angle( X_1(kmax-n:kmax+n)' * X(kmax-n:kmax+n) )

 

X_1 - предидущий выход FFT (вектор столбец),

X - текущий выход FFT,

kmax - индекс бина с максимальной энергией,

n - число 0...3

 

Можно бы и по всем. только не все бины одинаково полезны. Поэтому это будет не оптимально.

Усреднять надо в соответствии с уровнем сигнала каждого бина. Вычисление фазы в бинах с низким уровнем сигнала по отношению к шуму бесмысленны - это уже фаза компоненты шума а не сигнала. А уровень сигнала убывает при отклонении от центральных бинов достаточно быстро (быстрее чем 1/n), в то время как шум остается постоянным. Поэтому при низком уровне сигнал/шум лучше усреднять по 2-3-м, причем с учетом энергии (можно даже расписать точные формулы, выписывая плотность вероятности ошибки фазы в бине как функцию отношения сигнал/шум - это можно найти в учебниках по спектральному анализу или статрадиофизике - и максимизмизируя произведение вероятностей).

 

При очень высоком уровне сигнал/шум, хотя бы 20дб, действительно усреднять вроде можно по всем четырём (FFT_N/FFT_STEP в общем случае перекрытия, но тоже с учетом энергии бина), поскольку ошибки измерения фазы статистически независимы, но различны для разных бинов, в зависимости от энергии

 

да, f(k) - измеренная частота для бина k. Я не максимизировал функцию вероятности, а просто написал на вскидку формулу, которая будет вести себя оптимально для частот и посередине бина (f(k0)) и на границе бинов (f(k0)+f(k1))/2,

если этого не сделать, то точность упадёт на границе

 

ЗЫ. Вообще-то вычисление фазы через квадратуры всегда усиливает шум и помехи. Поэтому при вычислении частоты через фазу не стоит пренебрегать использованием спектральных окон, как в примере - чтобы подавить хвосты от возможных синусовых помех. Они и всегда мешают, а этому методу будут особенно

[alex_os 0]

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

 

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

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

[fontp 0]

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

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

 

Виноват. Про angle (в Матлабе?) не знал, думал среднее арифметическое фазы.

[alex_os 0]

Виноват. Про angle (в Матлабе?) не знал, думал среднее арифметическое фазы.

Сорри, не написал что матлаб. angle -аргумент комплексного числа.

 

[rudy_b 1]

Можно бы и по всем. только не все бины одинаково полезны...

Более правильно - по МНК подбирается параметры гаусса (огибающей) (при соответствующей весовой функции) с весами, пропорциональными мощности. При этом точность возрастает на порядок.