На ночь глядя наваял частотометр на C. Он слушает /dev/dsp, анализирует спектр и показывает самую сильную гармонику. Точность измерения здесь зависит от LENGTH, чем дольше замер, тем точнее определение. Чем ниже sample rate и короче длина куска, тем меньше жрёт памяти.

Я это сделал для замера сетевой частоты. В принципе, она должна быть 50 Hz, но иногда колеблется. Прямо сейчас я наблюдаю колебания в пределах 0.1 Hz при длине измерения 60 сек. Пока подключил микрофон, из-за этого вместо сетевого фона частенько ловится инфразвук от ветра между домами. Потом приделаю "измерительный трансформатор", показ времени и прикручу к gnuplot. Пусть по cron'у графики рисует и на сайте показывает.

Компилирую так:

cc -O2 50hz.c -lgsl -lcblas -lm

Работает под Linux, должно идти без переделки под FreeBSD и NetBSD. В данном варианте -- i386 only, потому что делает предположения о типах и порядке байтов.

/*
* 50hz.c
*/

#include <unistd.h>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/ioctl.h>
#include <stdlib.h>
#include <stdio.h>
#include <sys/soundcard.h>

#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>

#define LENGTH 60 /* how many seconds of hum to store */
#define RATE 1000 /* the sampling rate */
#define SIZE 16 /* sample size: 8 or 16 bits */
#define FMT AFMT_S16_LE /* change to AFMT_S16_BE for PPC and friends */

/* this buffer holds the digitized audio */
unsigned char buf[LENGTH * RATE * SIZE / 8];
signed short *sbuf;
double dbuf[LENGTH * RATE];

int main()
{
int fd; /* sound device file descriptor */
int arg; /* argument for ioctl calls */
int status; /* return status of system calls */
int i;
gsl_fft_real_wavetable * wavetable;
gsl_fft_real_workspace * workspace;
double freq, peak_freq;
double level, peak_level;

wavetable = gsl_fft_real_wavetable_alloc(RATE * LENGTH);
workspace = gsl_fft_real_workspace_alloc(RATE * LENGTH);
sbuf = (signed short *) buf;

/* open sound device */
fd = open("/dev/dsp", O_RDWR);
if (fd < 0) {
perror("open of /dev/dsp failed");
exit(1);
}

/* set sampling parameters */
arg = SIZE; /* sample size */
status = ioctl(fd, SOUND_PCM_WRITE_BITS, &arg);
if (status == -1) {
perror("SOUND_PCM_WRITE_BITS ioctl failed");
exit(1);
}
if (arg != SIZE) {
perror("unable to set sample size");
exit(1);
}

arg = 1; /* mono */
status = ioctl(fd, SOUND_PCM_WRITE_CHANNELS, &arg);
if (status == -1){
perror("SOUND_PCM_WRITE_CHANNELS ioctl failed");
exit(1);
}
if (arg != 1) {
perror("unable to set mono");
exit(1);
}

arg = RATE; /* sampling rate */
status = ioctl(fd, SOUND_PCM_WRITE_RATE, &arg);
if (status == -1){
perror("SOUND_PCM_WRITE_RATE ioctl failed");
exit(1);
}
if(arg != RATE){
perror("cannot set rate");
fprintf(stderr, "requested: %i set: %i ", RATE, arg);
exit(1);
}

arg = FMT; /* sample format */
status = ioctl(fd, SNDCTL_DSP_SETFMT, &arg);
if (arg != FMT){
perror("SNDCTL_DSP_SETFMT ioctl failed");
exit(1);
}

while (1) { /* loop until Control-C */
//printf("Say something: ");
status = read(fd, buf, sizeof(buf)); /* record some sound */
if (status != sizeof(buf)) {
perror("read wrong number of bytes");
exit(1);
}

/* convert 16-bit samples to double */
for (i = 0; i < LENGTH * RATE; i++) {
dbuf[i] = (double) sbuf[i];
//printf("%i %f ", i, dbuf[i]);
}

/* spectrum analyzer */
gsl_fft_real_transform(dbuf, 1, RATE * LENGTH, wavetable, workspace);
peak_freq = 0;
peak_level = 0;
for (i = 1; i < (LENGTH * RATE - 1); i+=2) {
freq = (i + 1.0) / (LENGTH * 2.0);
level = sqrt(dbuf[i]*dbuf[i] + dbuf[i+1]*dbuf[i+1]);
/* printf("%f %f ",
((double) i) / (LENGTH * 2.0),
sqrt(dbuf[i]*dbuf[i] + dbuf[i+1]*dbuf[i+1])
); */
if(level > peak_level){
peak_level = level;
peak_freq = freq;
}
}
printf("%f ", peak_freq);

//exit(0);
}
}

level = sqrt(dbuf[i]*dbuf[i] + dbuf[i+1]*dbuf[i+1]); это не совсем правильно, и я думаю Вы не инфразвук ловите, а результат именно этого...

p.s. Я не силен в С, но так уровень спектральных компонент не определяется особенно когда Вам нужно будет более менее реальную картину получить, а не условную....


p.p.s. А в прочем не берите в голову... х.з. может все и правильно... С не мой язык... Вам то всяко видней.

Этот кусок вычисляет модуль комплексного числа по реальной и мнимой части. Тут всё должно быть верно.

Вопрос в другом. Вот прибор:

http://www.power.kharkiv.com/ichps1.html

Как он измеряет с точностью 0.005 Гц за 1 сек? Каков принцип? Я тоже так хочу!

dimss

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

Как он измеряет с точностью 0.005 Гц за 1 сек? Каков принцип?
Да это же старо как мир..., по периоду сигнала... тут другая проблема точно и с нужным количеством знаков период измерить... хрен редьки не слаще...

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

1. Почему не даёт?
2. Что даёт?

dimss


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

> но если Вы захотите вычислить реальную ее величну допустм в дб, то такой подход даст ложный результат

А, это понятно. Но это и не нужно. Нужно только обнаружить сильнейшую гармонику.

Насчёт измерения периода -- мысль интересная. Грубо можно рассуждать так. Звуковушка может измерять с частотой 48000 измерений в секунду. Значит, один период может быть измерен за секунду с дискретностью до 1/960 его длительности. Сейчас за секунду достигается дискретность до 1/50...

dimss

Ну я бы посчитал вот так, период 50 Гц равен 0.02 сек, допустим измеряем с точностью +/- 0.001 тогда получаем точность измерения частоты 1/0.021 = 47.6 Гц это очень грубо, увеличим точность измерения периода до 0.0001 получаем 1/0.0201 = 49.75 точность примерно 0.3 герца, при этом измерить изменение длительности в 0.0001 сек частотой менее чем 1 Мгц не выйдет... вот собственно и все, это грубая и примитивная прикидка тем не менее указывает на сложность проблемы всего лишь с регистрацией периода. На самом деле проблем больше...

> при этом измерить изменение длительности в 0.0001 сек частотой менее чем 1 Мгц не выйдет

А почему, собственно? Звуковуха может давать отсчёт каждые 20 мксек, За секунду набегает около 50 периодов. Средний период вполне можно измерить с точностью примерно до десятков микросекунд. Конечно, нам будут мешать гармоники и помехи, но ведь сигнал можно отфильтровать, оставив только низкочастотные составляющие. И точность будет выше, чем у спектрального сособа.

Между тем набежали первые результаты утреннего измерения. Время одного измерения -- 120 сек. Сигнал снимался уже не с микрофона, а с куска провода.

50.050000
50.050000
50.041667
50.050000
50.066667
50.066667
50.058333
50.050000
50.041667
50.041667
50.041667
50.041667
50.050000
50.058333
50.050000
50.050000
50.050000
50.058333
50.058333
50.058333
50.058333
50.066667
50.050000
50.041667
50.050000
50.058333
50.050000
50.050000
50.041667
50.033333
50.033333
50.041667
50.041667
50.033333
50.041667
50.066667
50.058333
50.066667
50.058333
50.050000
50.041667
50.058333
50.050000
50.041667
50.041667
50.050000
50.058333
50.058333
50.050000
50.066667
50.050000
50.058333
50.058333
50.058333
50.066667
50.058333
50.058333
50.050000
50.041667
50.058333
50.058333
50.041667

А почему, собственно? В смысле? Хм... это железные законы бытия...

p.s. Ну да... я там малость обсчитался(с нулями запутался), для тех величин нужна частота 10 Кгц ну и соответственно на порядок увеличивать с каждым увеличением точности, от карты при 0.00002 cек точности получите точность измерения частоты примерно 50-(1/0.02002) ~ 0.05 Гц.

Сделал частотометр на принципе измерения периода. Быстродействие и точность повысились. О погрешности можно спорить, но при 5 сек. и 48000 kHz картина изменений частоты сети вполне ясна.

Принцип таков. Берём кусок звука и оставляем только полосу 45-65 Hz (ацкий FFT фильтр). Затем считаем фронты волн и замеряем время между первым и последним найденным фронтом. Вычисляем частоту. Повторяем.

Осталось приделать запись в файл и построение графиков для web.

---------------------------------

/*
* 50hz.c
*/

#include <unistd.h>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/ioctl.h>
#include <stdlib.h>
#include <stdio.h>
#include <sys/soundcard.h>

#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>

#define LENGTH 5 /* how many seconds of hum to store */
#define RATE 48000 /* the sampling rate */
#define SIZE 16 /* sample size: 8 or 16 bits */
#define FMT AFMT_S16_LE /* change to AFMT_S16_BE for PPC and friends */

/* this buffer holds the digitized audio */
unsigned char buf[LENGTH * RATE * SIZE / 8];
signed short *sbuf;
double dbuf[LENGTH * RATE];

int main()
{
int fd; /* sound device file descriptor */
int arg; /* argument for ioctl calls */
int status; /* return status of system calls */
int i;
gsl_fft_real_workspace * workspace;
gsl_fft_real_wavetable * r_wavetable;
gsl_fft_halfcomplex_wavetable * hc_wavetable;
double freq;
/* indices of 45 and 65 Hz for bandpass FFT filter */
int f1_index = 45 * LENGTH * 2 - 1;
int f2_index = 65 * LENGTH * 2 - 1;
int first_wave, last_wave, wavecount;

/* initialize fft engine */
workspace = gsl_fft_real_workspace_alloc(RATE * LENGTH);
r_wavetable = gsl_fft_real_wavetable_alloc(RATE * LENGTH);
hc_wavetable = gsl_fft_halfcomplex_wavetable_alloc(RATE * LENGTH);

sbuf = (signed short *) buf;

/* open sound device */
fd = open("/dev/dsp", O_RDWR);
if (fd < 0) {
perror("open of /dev/dsp failed");
exit(1);
}

/* set sampling parameters */
arg = SIZE; /* sample size */
status = ioctl(fd, SOUND_PCM_WRITE_BITS, &arg);
if (status == -1) {
perror("SOUND_PCM_WRITE_BITS ioctl failed");
exit(1);
}
if (arg != SIZE) {
perror("unable to set sample size");
exit(1);
}

arg = 1; /* mono */
status = ioctl(fd, SOUND_PCM_WRITE_CHANNELS, &arg);
if (status == -1){
perror("SOUND_PCM_WRITE_CHANNELS ioctl failed");
exit(1);
}
if (arg != 1) {
perror("unable to set mono");
exit(1);
}

arg = RATE; /* sampling rate */
status = ioctl(fd, SOUND_PCM_WRITE_RATE, &arg);
if (status == -1){
perror("SOUND_PCM_WRITE_RATE ioctl failed");
exit(1);
}
if(arg != RATE){
perror("cannot set rate");
fprintf(stderr, "requested: %i set: %i ", RATE, arg);
exit(1);
}

arg = FMT; /* sample format */
status = ioctl(fd, SNDCTL_DSP_SETFMT, &arg);
if (arg != FMT){
perror("SNDCTL_DSP_SETFMT ioctl failed");
exit(1);
}

while (1) {
status = read(fd, buf, sizeof(buf));
if (status != sizeof(buf)) {
perror("read wrong number of bytes");
exit(1);
}

/* convert 16-bit samples to double */
for (i = 0; i < LENGTH * RATE; i++)
dbuf[i] = (double) sbuf[i];

/* bandpass FFT filter */
gsl_fft_real_transform(
dbuf, 1, RATE * LENGTH, r_wavetable, workspace);
for (i = 0; i < (LENGTH * RATE); i++)
if((i < f1_index) || (i > f2_index)) dbuf[i] = 0;
gsl_fft_halfcomplex_transform(
dbuf, 1, RATE * LENGTH, hc_wavetable, workspace);

/* let's count waves! */
wavecount = 0;
first_wave = 0;
last_wave = 0;
for (i = 1; i < LENGTH * RATE; i++) {
if((dbuf[i-1] < 0) && (dbuf[i] >= 0)){
if(!wavecount) first_wave = i;
last_wave = i;
wavecount++;
}
}

/* finally calculate mains frequency */
freq = ((double) RATE) / ((last_wave - first_wave) / (wavecount - 1.0));
printf("%f ", freq);
}
}


--------------------------------------------

Компилируем:

cc -Wall -O2 -s -o 50hz 50hz.c -lgsl -lgslcblas -lm

-----------------------------------------------

Примерный вывод:

49.995399
49.995608
49.995817
49.996654
49.997072
49.997490
49.997908
49.999372
49.998954
49.997490
49.997908

Это изменение частоты в офисе за последние полчаса.

http://www.solutions.lv/dimss/01.png

Осталось сделать автоматическое рисование графиков с часами и минутами.

Ну вот, поставил на автомат:

http://www.solutions.lv/50hz/

А что за библиотека для БПФ?

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

Mixa
А что за библиотека для БПФ?

GSL -- GNU Scientific Library. Удобнейшая штука!

Valery
А где гарантия в стабильности кварцевого генератора звуковухи, от которой зависит точность задания перида счета импульсов?

Нигде. Чем богаты, тем и рады. Но, как говорится, для наших целей -- достаточно и Creative SB PCI 128. При необходимости можно поставить и хороший измерительный ADC вместо звуковухи. Мне интересно было саму программу написать.

Температурная нестабильность частоты генератора, нестабильность питающего напряжения и т.д. звуковухи могут дать изменение частоты в 0,001 процента (это реальная величина), а измеренные Вами отклонения частоты (в процентах) укладываются в этот диапазон. Вот найти бы эталон частоты 50 Гц!
Как вариант: настраиваем приемник на эталонную частоту в режиме SSB c разносом 100- 200 Гц и проводим калибровочные измерения.

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

Код нынешнего варианта можно получить здесь:

http://dimss.solutions.lv/soft/50hz.tgz

Тем временем в моём воображении зреет следующая стадия проекта -- объединение первого и второго метода. Тогда можно будет измерять частоту сигнала, лежащего под помехами, если только спектр помех более размазан.

Внесены изменения в программный частотометр для Linux:

* автоматическая настройка полосового фильтра на основную частоту входного сигнала;
* вывод в ежедневный файл теперь делается отдельной функцией;
* компилируется без замечаний GCC 4.0.1.

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

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

К сожалению, сейчас мне негде проводить непрерывные наблюдения за частотой электросети. Вся надежда на добровольцев, располагающих постоянно работающим Linux-компом со звуковухой и Web-сервером.

Код частотометра можно взять здесь:

http://dimss.solutions.lv/soft/50hz.tgz

На выходных я взял свой старый код и переделал само измерение частоты. Теперь частота измеряется при помощи алгоритма ФАПЧ, т.е. воспроизводится частота около 50 Гц и сравнивается по фазе со входным сигналом. Результат сравнения интегрируется и используется для управоления воспроизводимой частотой. Когда происходит захват, воспроизводимая частота является оценкой частоты входного сигнала.

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

Фильтр я пока оставил старый (ДПФ), он не ухудшает точность при новом методе. Но в дальнейшем думаю перейти на более простой БИХ-фильтр, режущий все выше 100 Гц и убирающий долгосрочную постоянную составляющую. Думаю, этого будет достаточно при наличии измерительного трансформатора (220 на 0.5 в линейный вход).

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

http://stalin.lv/soft/50hz-20080102.tgz

Запустите, пожалуйста, по крону - интересно графики смотреть :)

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

Основные нововведения:

1) добавлен второй канал сравнения фазы, в котором воспроизводимый сигнал сдвинут на 90 градусов относительно первого канала, т.е. синус-косинус. Если в первом канале при захвате частоты на выходе фильтра устанавливается величина около 0.5, то во втором она должна отличаться точно на 0.5 от первой. Причем отличие от идеального значения должно быть у обоих каналов одинаковое. Следовательно, при точном удержании фазы "разница разниц" фаз будет стремиться к нулю. Это можно использовать для определения правдоподобия выходной величины в данный момент.
2) Измеряется не частота непосредственно, а ее 7-я гармоника. Для этого сигнал фильтруем 45-55 Гц, делаем прямоугольным, снова фильтруем 315-385 Гц. Измеряем. Результат делим на 7.

Так у нас повышается точность и предотвращается выдача некорректных результатов до полного захвата (т.е. во время колебаний в цепи ФАПЧ).

Завтра выложу свежие исходнички. Руки все не доходят сделать измерительный трансформатор, чтобы действительно поставить измерительный процесс на недельку. С графиками.