[ivan219 0]

1. Обнаружена методическая погрешность из-за использования оконной функции. Погрешность имеет синусоидальный характер с периодом пол-бина

Не совсем синусоида смотрите рисунок. Причём чем больше оконная функция имеет подавление второго лепестка, тем точнее результат, это справедливо для чистой синусоиды. Но для зашумленного сигнала нужно подбирать оконную функцию по максимальным значениям. Ещё можно сделать таблицу погрешности и тем самым улучшить результат. Так же было замечено, что результат всегда находится ближе к центральному бину чем реальное значение. И можно ещё более точно определить значение, если умножить на комплексную экспоненту и снести ближе к центральному бину.

Изменено 25 марта, 2011 пользователем ivan219

[Ruslan1 16]

Не совсем синусоида смотрите рисунок. Причём чем больше оконная функция имеет подавление второго лепестка, тем точнее результат, это справедливо для чистой синусоиды. Но для зашумленного сигнала нужно подбирать оконную функцию по максимальным значениям. Ещё можно сделать таблицу погрешности и тем самым улучшить результат. Так же было замечено, что результат всегда находится ближе к центральному бину чем реальное значение. И можно ещё более точно определить значение, если умножить на комплексную экспоненту и снести ближе к центральному бину.

Ага, согласен. и у меня то же самое, даже график приводить не буду- один в один. конечно не синусоида, под "синусоидальностью" я подразумевал форму, издалека похожую на синус :) Ну и период равен бину, а не пол-бина, как в моем сообщении, это я куда-то не туда посмотрел.

 

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

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

 

Вот привожу файлик исходных семплов АЦП (странно, файл с расширением csv загружать запрещено, пришлось в txt переобозвать).

ADC.TXT

 

Файл имеет простой текстовый csv-формат, который хоть в ексель хоть в матлаб подсовывается с пол-пинка. на всякий случай привожу код, как его в матлаб прочитать:

Fs = 31250
FFT_N_SAMPLES = 4096
%% input data
% read data from file
filename = 'ADC.TXT';
fid=fopen(filename','rt'); 
data_in = fscanf(fid, '%d,');
fclose(fid);

 

Мои результаты вычислений:

1. максимальный бин (индексы начинаются с единицы, как принято в Матлабе):

n_bin_max = 109

 

2. расстояние от максимального бина:

d = -0.3976

 

3. собственно координаты в бинах :

bin = n_bin_max + d = 108.6024

 

4. Частота при указанных выше 31250 SPS и 4096 точках:

freq2 = (bin-1) * Fs / (FFT_N_SAMPLES) = 820.9413

 

Согласно матлабу, результат может изменится на величину до -0.02 Гц при изменении коэффициента альфа в Гауссовском окне (от 4 до 10). Но в любом случае не достигается величина, показываемая методом подсчета пересечений через ноль (разница до -0.09 Гц от FFT-based метода)

 

Кстати, у меня так и получается, что "результат всегда находится ближе к центральному бину чем реальное значение". Насчет "всегда" не знаю, но вот та же петрушка.....

 

[Alex11 3]

Поанализировал я тут Ваш файлик. Очень тяжело искать черную кошку в темной комнате, особенно, еслиее там нет. Так вот. Мало того, что это не синусоида, мало того, что зашумлено, но частота сигнала изменяется на протяжении этого семпла. Так о чем может идти речь? На 1/4 семпла от начала частота 820.9137, на 1/2 - 820.9203, на 3/4 - 820.9486. За последний знак не поручусь, т.к. частота плывет непрерывно и это тоже усредненные данные. Интегрально мой алгоритм дает частоту 820.9370. И это естественно, т.к. окно отрезает начало и конец, затем Фурье дает некоторое интегральное значение.

А откуда значение 820.9413, которое хотелось получить?

[Ruslan1 16]

Поанализировал я тут Ваш файлик. Очень тяжело искать черную кошку в темной комнате, особенно, еслиее там нет. Так вот. Мало того, что это не синусоида, мало того, что зашумлено, но частота сигнала изменяется на протяжении этого семпла. Так о чем может идти речь? На 1/4 семпла от начала частота 820.9137, на 1/2 - 820.9203, на 3/4 - 820.9486. За последний знак не поручусь, т.к. частота плывет непрерывно и это тоже усредненные данные. Интегрально мой алгоритм дает частоту 820.9370. И это естественно, т.к. окно отрезает начало и конец, затем Фурье дает некоторое интегральное значение.

А откуда значение 820.9413, которое хотелось получить?

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

 

Насчет несинусоидальности. Для меня струнные датчики новая область, но везде пишут что это именно синус на выходе. И датчик у меня честный промышленный, не самоделка. значит, что-то с условиями я не соблюдаю у себя на столе. Ну или как вариант- моделируется это колебание чем-то(дом шатается, я уже такое проходил с сейсмосенсорами, 7-й этаж уже болтается вполне предсказуемо), а датчик чувствительный. Это же относительно изменения частоты- может это так и должно быть, реально меняется показание.

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

Но в данный момент очень хочется попробовать доползти до Вашей цифры.

 

Как я уже писал, для меня ЦОС это новая область. Я уже не говорю о том, что Матлаб впервые установил дней десять назад.

До всего дохожу в рабочем порядке, читая форумы и интернет-ресурсы. В результате родился этот файл testing_on_real_data2.txt (нужно переименовать в *.m). ниже привожу его же.

 

он и возвращает указанное значение: Frequency_Response = 820.9413

 

FFT_N_SAMPLES = 4096;
Fs =31250;
alpha = 3.9
    
%% input data
% read data from file
filename = 'ADC.CSV';
fid=fopen(filename','rt'); 
data_in = fscanf(fid, '%d,');
fclose(fid); 
%% windowing
for i=1:FFT_N_SAMPLES
        %Gaussian:
        tmp_n = abs(i - (FFT_N_SAMPLES)/2);
        tmp = alpha * tmp_n / (FFT_N_SAMPLES/2);
        coeff = exp(-0.5 * tmp*tmp);
        data_win(i) = data_in(i) * coeff;
end
%% fft
data_fft = fft (data_win,FFT_N_SAMPLES);
%% magnitude
for i=1:FFT_N_SAMPLES/2
    data_fft_abs(i) = 2.0 * abs(data_fft(i));
end
%% maximum searching
max = data_fft_abs(1);
n_bin_max = 1;
for i=1:FFT_N_SAMPLES/2
    if (max < data_fft_abs(i))
        max = data_fft_abs(i);
        n_bin_max = i;
    end
end
%% Gaussian interpolation
            m1 = data_fft_abs(n_bin_max-1);
            m2 = data_fft_abs(n_bin_max);
            m3 = data_fft_abs(n_bin_max+1);
            top_part = log (m3/m1);
            bot_part = 2*log(m2*m2/(m3*m1));
            d = top_part/bot_part;
% f = ( i + d ) * ( sr / n )
            n_bin_response = n_bin_max + d;
            Frequency_Response = (n_bin_response-1) * Fs / (FFT_N_SAMPLES)

 

Самым слабым местом мне кажется оконная функция, хотя может еще где напортачил.

 

Поанализировал я тут Ваш файлик. .... Мало того, что это не синусоида, мало того, что зашумлено..........

Сейчас нашел файлик, который я оцифровывал через саундкарту с этого же датчика в тех же условиях, но с другим аналоговым трактом и другим АЦП. Все-таки больше похоже на синус. in_data32000SPS_4096pnt.zip

Файл читается в Матлаб после распаковки командой

[data,time] = daqread('in_data32000SPS_4096pnt.daq');

Значит я не только в математике напортачил..... Это я конечно завтра начну рыть....

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

[Alex11 3]

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

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

[ivan219 0]

Поможет БИХ фильтр. Если есть вычислительные ресурсы, то хороший фильтр избавит от гармоник. Единственно, что плохо так это вычисление в несколько этапов.

[Alex11 3]

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

Теперь про результаты анализа файла с 32к оцифровки. Тяжело плавать в серной кислоте с оторванными ногами. У меня матлаб не стоит, так что пришлось потрошить файл вручную. Где-то я ошибся на один семпл и у меня прошел сбой по частоте. Но в целом этот файл гораздо лучше последнего, частота болтается не больше, чем 2 единицы второго знака после запятой, хотя амплитуда маленькая и шумов довольно много. Если сохраните файл в текстовом виде, могу посчитать еще раз поточнее. Первая прикидка частоты - 822.557.

[Ruslan1 16]

Теперь про результаты анализа файла с 32к оцифровки. Тяжело плавать в серной кислоте с оторванными ногами. У меня матлаб не стоит, так что пришлось потрошить файл вручную. Где-то я ошибся на один семпл и у меня прошел сбой по частоте. Но в целом этот файл гораздо лучше последнего, частота болтается не больше, чем 2 единицы второго знака после запятой, хотя амплитуда маленькая и шумов довольно много. Если сохраните файл в текстовом виде, могу посчитать еще раз поточнее. Первая прикидка частоты - 822.557.

Да конечно пожалуйста. Сконвертировал из *.daq и добавил в архив и текстовый файл, расширение "csv". in_data32000SPS_4096pnt.zip

 

И опять засада: у меня 822.5698 Hz :(

 

четверти от выборки (по 1024 точек):

1: 822.5570

2: 822.5868

3: 822.5385

4: 822.5531

половинки от выборки (по 2048 точек):

1: 822.5842

2: 822.5649

Полная выборка (4096 точек): 822.5698

 

Получается, что измеряемая частота сильно плывет? такое теоретически возможно. Но как-то не вяжется с тем что я вижу глазами: я последовательно применяю два метода (подсчет периодов методом пересечение нуля и FFT-based: первый-второй-первый-второй...), и вычисленная каждым из них величина колеблется около своего значения частоты, причем FFT-baised больше. И разность того же порядка как разница между Вашей и моей цифрами. :(

[Alex11 3]

Вот уточненные результаты: средняя частота 822.5666. Если смотреть по частям (я, правда, делал несколько иначе - оставлял FFT на 4096, но делал узкое окно и сдвигал его относительно данных), то на 1/4 от начала частота 822.5897, на 1/2 - 822.5614, на 3/4 - 822.5633. При делении на 8 получается следующее:

822.5477

822.6066

822.5791

822.5578

822.5175

822.5213

822.5499

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

Я посмотрел на модели - синус плюс белый шум примерно равный тому, что есть в этом файле. Ошибка частоты - примерно 0.03 Гц максимум при делении на 8. Таким образом, здесь кроме шума действительно есть небольшое плавание частоты (примерно на те же 0.03 Гц). На полном интервале данный шум приводит к ошибке меньше 0.01 Гц. Отсюда следуют рекомендации (очевидные) - увеличивайте отношение сигнал/шум и боритесь с плаванием частоты.

 

[Ruslan1 16]

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

Я посмотрел на модели - синус плюс белый шум примерно равный тому, что есть в этом файле. Ошибка частоты - примерно 0.03 Гц максимум при делении на 8. Таким образом, здесь кроме шума действительно есть небольшое плавание частоты (примерно на те же 0.03 Гц). На полном интервале данный шум приводит к ошибке меньше 0.01 Гц. Отсюда следуют рекомендации (очевидные) - увеличивайте отношение сигнал/шум и боритесь с плаванием частоты.

Понял, большое спасибо!

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

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

 

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

[Alex11 3]

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

[Ruslan1 16]

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

 

на чистом синусе все шикарно, внутри бина ошибка (разность вычисленной и теоретической частот) соответствует предсказанной для гауссовского окна, именно так я и подбирал коэффициенты для окна (чтобы не более 0.001 Гц максимальная ошибка была).

Более того- и в приборе при подаче чистого синуса все отлично, ну прямо почти единицы наносекунд совпадают с частотомером.

А в реальности что-то не очень. Одни только очень серьезные палки на нечетных гармониках чего стоят. А фильтрануть не могу- полоса прибора позволяет куче гармоник пролазить. Реально хорошо вижу все нечетные гармоники, сколько аппаратный ФНЧ оставляет.

 

С шумами не смотрел на модели вообще, есть такой грех. Смотрел только разрешающую способность и точность при чистой синусоиде и радовался как слон :) Посмотрю... Кстати и с гармониками посмотрю....

 

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

Кстати я приводил картинку, которую производители суперкрутого прибора привели:вот. Ну нет у них гармоник, и сигнал затухает в сотни раз быстрее чем у меня. Может еще и спецдатчики нужны для такого прибора.... Или демпфирование.....

 

Но это все так, рюшечки для себя любимого и набирание опыта с Матлабом. Поставленная при старте топика задача решена :)

 

Еще раз спасибо!

:)

[Alex11 3]

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

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

[Ruslan1 16]

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

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

Попробовал демпфировать датчик резистором- получил увеличение отношения сигнал/гармоника в 2 раза! Это я катушку датчика конкретно так нагружаю.... бум думать, тоже интересно...... Но не работает никто в токовом режиме, это я уже патент какой могу заявлять наверное :)

Гармоники проверил (на модели). Действительно зря боялся. На определение основной частоты при выбранном мной окне не влияют совсем: при 3-й,5-й 7-й гармониках с амплитудой равной основной гармонике рассчетная частота отличается от поданной примерно на 0.0006 Гц. То есть с этой стороны все в порядке.

 

А насчет той красивой картинки- так они точность 0.013% от значения обещают. даже для 100Гц это 0.013Гц, для более реальных 500Гц и выше- 0.065 Гц :) Не лучше чем у меня :)

[Alex65111 0]

Долго читал про потенциальную точность измерения, про то что здесь получается у народа (доли герца), но что-то с моими данными не могу сообразить.

 

 

1 вопрос. Имеются данные dat1.zip (после распаковки в матлабе загружаются load dat1). Как для данного конкретного случая (сигнал/шум, объем выборки, частота дискретизации 68кГц) оценить потенциальную точность измерения частоты? Формулу в предыдущих постах видел, но не пойму какое сигнал/шум подставлять, то ли во всей полосе, то ли только в полосе сигнала.

 

2. Насколько оптимально по точности измерение частоты следующим подходом, какую точность здесь можно ожидать?

 

Fs=68000;

N=512;

df=Fs/N;

 

sp=fft(s,N);

 

fmax=find(abs(sp(1:length(sp)/2))==max(abs(sp(1:length(sp)/2))));

 

f_ocen=(fmax(1)-1)*df;

 

f_s=f_ocen-df/1;

f_e=f_ocen+df/1;

 

m = 128*2;

w = exp(j*2*pi*((f_s-0)-(f_e+0))/(m*Fs));

a = exp(j*2*pi*(f_s-0)/Fs);

 

L1 = czt(s,m,w,a);

f_izm=find(L1==max(L1))*(2*df/m)+f_s;