как это сделано: FFT4096, 100-6500Hz, freq.resolution 0.001Hz, как достигнуто такое разрешение по частоте? Опции

[Ruslan1]

Здравствуйте!

Помогите пожалуйста, спать не могу, все думаю......

Накопал в интернете, не понимаю как сделано. а хочется понять и применить методу.

короткая листовка: http://www.campbellsci.com/documents/produ...es/b_avw200.pdf
полный мануал: http://www.campbellsci.com/documents/manuals/avw200.pdf


Прибор определяет частоту синусоидального сигнала, который может быть в диапазоне 100 - 6500 Гц, с разрешением 0.001Гц и точностью 0.013% от измеренной величины. Указано что используется 4096-точечное FFT. Интересно, что еще кроме упомянутого FFT применяется для получения такого разрешения по частоте?

Как это может быть сделано? Взяли FFT, определили частоту грубо, потом детально обнюхали область вокруг найденной частоты, получили частоту точно?


Что абсолютно точно известно (это физика процесса):
1. Полезный сигнал- честная одночастотная синусоида (дернули струну, потом ищут ее частоту собственных колебаний. то есть это затухающее колебание механической струны, которое преобразовано в электрический сигнал). Шумы возможны.
2. Общая длительность исследуемого сигнала доли секунды (ну пусть 300 ms)


Что говорит изготовитель (интересный текст и какие-то картинки начинаются со страницы 85, еще есть спецификация на странице 9):
1. The AVW200 uses an audio A/D for capturing the sensor’s signal. The number of samples acquired in this period is 4096 points. A Fast Fourier Transform (FFT) algorithm is used to create a frequency spectrum.
2. the spectral analysis gives improved frequency resolution (0.001 Hz rms) during quiet conditions.
3. the AVW200 Fourier transforms the recorded response and analyzes the resulting spectrum to determine the wire’s resonant frequency. This analysis also provides diagnostic information indicating the quality of the resonant-frequency measurement.
4. Measurement resolution: 0.001 (Hz RMS)
5. Accuracy basic: +/- 0.013% of reading.
6. DF measurement time between 1.6 to 1.8 second
7. SYSTEM: PROCESSOR: Hitachi H8S 2324 (16-bit CPU with 32 bit internal core), MEMORY: Either 128 or 512 kbytes of SRAM; 2 Mbyte of OS Flash
8. частоту проца не знаю, но пишут про дополнительный ток потребления во время измерения 27mA@12V. Судя по даташиту проца получается что используется он по максимуму, 25 МГц тактовой.


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

[ivan219]

И всё таки не много не понятно.
Что делает эта функция?

Код

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

[tmtlib]

И всё таки не много не понятно.
Что делает эта функция?

Код

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]

Фазу приводят к диапазону -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 элемента считаем фазу. Если по простому, то фаза[i]:=ArcTan(Im[i] / Re[i]); Можно сделать получше, расписав разные случаи (если Re[i]=0 и т.д.)
3) храним фазы за предыдущий проход FFT и за текущий. m_fftLastPhase[i] - это старая фаза, а phase - это на самом деле phase[i], т.е. фаза[i] от текущего 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[i] - это старая фаза, а 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]

...
(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]

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

CODE

        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,
если этого не сделать, то точность упадёт на границе

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

[fontp]

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

Так предложили уже взвешивать, если говорить о алгоритме с разностью фаз ДПФ перекрывающихся фреймов.
Или Вы о чём то своём?

Сообщение отредактировал fontp - Nov 16 2011, 11:10