Jump to content

    
Sign in to follow this  
Lost_Viking

Модифицированный алгоритм Герцеля

Recommended Posts

Взялся за реализацию. Прочитал про алгоритм тут: dsplib

Написал такую реализацию:

#define Fd 9000 //частота дискретизации
#define f1 600  //частота тона
#define f_search 600 //искомая алгоритмом частота
#define N 70 //размер массива
int main()
{
unsigned char ADC[350]; //Массив, эмулирующий поток из АЦП
unsigned char buff[N]; //Буфер, в который пишутся данные из АЦП
float I[2]; 
float Q[2];
float a;
float b;
float c = 2* M_PI *f1/Fd; //Заранее рассчитаем то, что не будет меняться
float module; //здесь будет лежать модуль комплексного числа
I[0]=0; I[1]=0; //инициализируем массивы
Q[0]=0; Q[1]=0;
for (int i =0; i<350; i++)
{
    ADC[i] = (unsigned char) 100*sin(c*i)+100; //заполняем "поток" АЦП синусом с амплитудой от 0 до 100

}

a=cos(2*M_PI*f_search/N); //рассчитываем к-ты
b=sin(2*M_PI*f_search/N);

for (int n=0; n<70; n++) //начальное заполнение буфера
{
    I[0]=I[1]; //Готовим предыдущее I
    Q[0]=Q[1]; 
    buff[n]=ADC[n]; // Просто заполняем буфер данными из АЦП
    I[1]=a*(I[0]+buff[n])-b*Q[0]; //Все по формуле 16 из статьи
    Q[1]=a*Q[0]+b*(I[0]+buff[n]);

}
for (int n=70; n<350; n++) //Когда буфер заполнился, то запускаем цикл с выталкиванием
{
    I[0]=I[1];
    Q[0]=Q[1];
    I[1]=a*(I[0]+ADC[n]-buff[0])-b*Q[0]; //В качестве i(n) выступает ADC(n), где n=70. Следовательно buff(n-N)=buff(0)
    Q[1]=a*Q[0]+b*(I[0]+ADC[n]-buff[0]); 
    for (int nn=0; nn<69; nn++) //перемещаем элементы массива на 1 влево
    {
        buff[nn]=buff[nn+1]; 
    }
    buff[69]=ADC[n]; //Заносим полученное значение АЦП на данном цикле работы алгоритма в последний эллемент буфера
    module = sqrt(I[1]*I[1]+Q[1]*Q[1]); //вычисляем модуль комплексного числа
}
//конец работы алгоритма
}

Для удобства чтения скрин из саблайма: screenshot

 

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

 

Скрин вывода программы переменной module:screenshot

 

Что я делаю не так?

Edited by Lost_Viking

Share this post


Link to post
Share on other sites

На елку такие скриншоты можно вешать, и гирлянды не надо. Все люди как люди, выкладывают код между соответствующими тэгами.

Share this post


Link to post
Share on other sites
На елку такие скриншоты можно вешать, и гирлянды не надо. Все люди как люди, выкладывают код между соответствующими тэгами.

Ок,понял. Через 30 минут исправлю

 

Share this post


Link to post
Share on other sites

В книге Р. Лайонса "Цифровая обработка сигналов" алгоритм Герцеля описан. Кажется, если сделать по формулам, то и работать должно безупречно.

Share this post


Link to post
Share on other sites

Да, он там описан. И тот алгоритм работает как надо. Но дело в том, что тот алгоритм работает после N семплов. То есть заполняем буфер, а потом обрабатываем его.

 

А на dsplib предложен алгоритм, который позволяет на лету обрабатывать семплы. Так вот разобраться с ним не могу.

 

Немножко напишу про свое понимание алгоритма, который на dsplib. Если перейти по ссылке сюда: dsplib, то можно увидеть вот что:

формула 12: goertzelmod_html_m2caca4a6.gif. О ней говорится, что это входной сигнал. На форуме dsplib мне ответили, что q(n) нужно принимать равными нулю. Вполне логично. Соответственно сама рекурсивная формула имеет вид:

goertzelmod_html_m7d342902.gif

Из нее убираем все q. Но тогда какой смысл вываливать такую формулу, если ее можно сократить? Это смущает

Edited by Lost_Viking

Share this post


Link to post
Share on other sites

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

 

Share this post


Link to post
Share on other sites
Наверное тогда лучше так: Скользящее ДПФ

 

Что-то я не понял каким образом здесь решается проблема конечной разрядной сетки, всё тот же компелексный множитель в интеграторе остаётся?

Share this post


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Sign in to follow this