Elektronik / Elektronik Kaynakları/

Teoriden Gerçeğe: Goertzel Algoritması (Vol. II)

Sponsorlu Bağlantılar

Daha önceki yazımda teorisinden bahsedip, bilgisayar ortamında örneklerini yaptığım Goertzel Algoritmasını bugün en sonunda gerçeğe dökebildim.

Uygulamada temel bileşen ile (50Hz) özellikle elektrik ve güç elektroniği ile ilgilenenlerin başını sıkça ağrıtan ilk 10 harmoniği (150-250-350-450-550-650-750-850-950-1050) algılamaya çalıştım. Yaptığım uygulama ayrıca PIC24F serileri ile örnek alımını ve bu serilerin ADC ve Timer birimlerini öğrenmek isteyenlere güzel bir örnek teşkil edeceğini düşünüyorum.

Yaptığım uygulamaya ait blok diyagram aşağıda görülebilir.

GOERTZEL

Yukarıdaki blok diyagramın görevini yerine getiren kodu aşağıda görebilirsiniz.

/*
 * Yazar: Fırat Deveci (FxDev)
 * Konu : Goertzel Algoritması
 */

#include "p24FJ64GA002.h"
#include "math.h"

_CONFIG1(JTAGEN_OFF & GCP_OFF & GWRP_OFF & BKBUG_OFF & COE_OFF & ICS_PGx1 & FWDTEN_OFF & WINDIS_OFF);
_CONFIG2(IESO_OFF & SOSCSEL_LPSOSC & WUTSEL_FST & FNOSC_FRCPLL & FCKSM_CSDCMD & OSCIOFNC_OFF & IOL1WAY_OFF & I2C1SEL_PRI & POSCMOD_NONE)

#define Fs              5000            // Örnekleme frekansımız, izleyeceğimiz en yüksek frekans değerinin en az iki katı
#define BINS            11              // 11 temel frekans izlenecek
#define Pi              3.141592654     // Pi sayısı tanımlanıyor
#define GOERTZEL_N      100             // 128 örnek alınacak
#define Level_Segment   7
#define divider         16000000/Fs

const int   freqs    [BINS] = {50, 150, 250, 350, 450, 550, 650, 750, 850, 950, 1050};
float       samples  [GOERTZEL_N];

float       coeffs   [BINS];
float       prev1    [BINS];
float       prev2    [BINS];
float       magnitude[BINS];
int         levels   [BINS];
float       n        [GOERTZEL_N];

unsigned char sample_counter    = 0;
unsigned char sample_complete   = 0;

void timer2_init(void)
{
    // 16Mhz/ftimer=PR2 olmalı.
    PR2                 = divider;
    T2CON               = 0x0000;     // 1:1, Fosc/2
    IPC1bits.T2IP       = 6;          // Kesme öncelik sırası 6
    IFS0bits.T2IF       = 0;          // Timer1 kesme bayrağı temizleniyor
    IEC0bits.T2IE       = 1;          // Timer1 kesmesi aktif
    T2CONbits.TON       = 1;          // Timer1 açılıyor
}

void ADC_init(void)
{
    AD1PCFG             = 0xFFFC;       // AN0 ve AN1 ADC portu, diğer AN portları digital I/O
    AD1CON2             = 0x0000;       // AVdd ve AVss seçildi, 28pinli modeller için bunlar zaten VDD ve VSS'ye bağlı
    AD1CON3bits.ADRC    = 0;            // Tcy kullanılacak
    AD1CON3bits.SAMC    = 3;            // Minimum Tad süresi 75ns'den büyük seçilmeli, 32MHz kristal için Tcy=62.5ns, bu değer için 2 yeterlidir
    AD1CON3bits.ADCS    = 2;            // ADCS=(Tad/Tcy)-1 şeklinde seçilir
    AD1CON1             = 0x00E0;       // Atomatik örnekleme ve çevrim yapılacak, ASAM yada SAMP bitine emirle çevrim başlayacak
    AD1CSSL             = 0x0000;       // Auto tarama yapılmayacak
    AD1CON1bits.ADON    = 1;            // ADC açılıyor
}

int ADC_read(unsigned char channel)
{
    AD1CHS              = channel;
    AD1CON1bits.DONE    = 0;
    AD1CON1bits.SAMP    = 1;           // Çevrim başlama komutu veriliyor
    while(AD1CON1bits.SAMP);
    while(!AD1CON1bits.DONE);          // Çevrimin bitmesi bekleniyor

    return ADC1BUF0;
}

// Goertzel için katsayılar hesaplanıyor

void coefficient(void)
{
    unsigned int k, i;
    for(i=0;i < BINS;i++)
    {
        k           =(unsigned int)(0.5 + (float)GOERTZEL_N * freqs[i] / Fs);
        coeffs[i]   = 2.0 * cos(2.0 * Pi * ((float)k/GOERTZEL_N));
    }
}

int goertzel(void)
{
    int big =0,
        i   =0;
    float val = 0,
          max = 0;

    // Geçici dizilerin içerisi temizleniyor
    for(i=0;i < BINS;i++)
    {
        prev2[i]=0;
        prev1[i]=0;
    }

    // GOERTZEL Algoritması İşletiliyor
    for(sample_counter=0;sample_counter < GOERTZEL_N; sample_counter++)
    {
        for(i=0;i < BINS;i++)
        {
            val     = coeffs[i] * prev1[i] - prev2[i] + (float)samples[sample_counter];
            prev2[i]= prev1[i];
            prev1[i]= val;
        }
    }

    // Genlikler belirleniyor
    for(i=0;i < BINS;i++)
    {
        magnitude[i] = (prev1[i] * prev1[i]) + (prev2[i] * prev2[i]) - (coeffs[i] * prev1[i] * prev2[i]);
        if(magnitude[i]>max)
        {
            max=magnitude[i];
            big=i;
        }
    }

    for(i=0;i < BINS;i++)
    {
        levels[i]=(int)(Level_Segment*(float)(magnitude[i]/max));
    }

    return big;
}

int main(void)
{
    CLKDIV=0;
    TRISB=0x0000;

    lcd_init();

    coefficient();
    ADC_init();
    timer2_init();

    lcd_gotoxy(2,1);
    lcd_yaz("FxDev.org");

    for(;;)
    {
        if(sample_complete)
        {
            goertzel();     // Goertzel hesaplanıyor
            lcd_gotoxy(1,1);
            veri_yolla(levels[0]);
            veri_yolla(levels[1]);
            veri_yolla(levels[2]);
            veri_yolla(levels[3]);
            veri_yolla(levels[4]);
            veri_yolla(levels[5]);
            veri_yolla(levels[6]);
            veri_yolla(levels[7]);
            veri_yolla(levels[8]);
            veri_yolla(levels[9]);
            veri_yolla(levels[10]);
            sample_complete = 0;
            sample_counter  = 0;
            IEC0bits.T2IE   = 1;
            T2CONbits.TON   = 1;
        }

    }
}

void __attribute__((__interrupt__, auto_psv)) _T2Interrupt(void)
{
    samples[sample_counter++]=(ADC_read(0)-512);
    if(sample_counter==GOERTZEL_N)
    {
        IEC0bits.T2IE   = 0;
        T2CONbits.TON   = 0;
        sample_complete = 1;
    }
    IFS0bits.T2IF = 0;
}

Yukarıdaki kodun daha önce burada anlattığım C kodlarıyla ne kadar fazla benzerlik gösterdiğini görebilirsiniz.
Uygulamaya ait çalışma videosunu aşağıda seyredebilirsiniz.

Bir başka çalışmada görüşmek üzere.
Herkese kolay gelsin.

  • murat

    Oldukça basit ve güzel bir uygulama olmuş. Bir çok farklı amaç için kullanılabilecek ve MCU kaynaklarını oldukça verimli kullanan bir algoritma.

  • Ali

    Güzel bir çalışma olmuş, tebrikler.

    Koddaki;
    >> #define divider 16000000/Fs
    ve
    >> samples[sample_counter++]=(ADC_read(0)-512);
    kısımlarını anlayamadım.
    Öncelikle MCU 16Mhz de çalıştığı için mi bölme değeri 16M oldu?
    Ve diğerinde neden adc değerinden 512 çıkarılıyor?

    Ayrıca adc okumasını kesme içinde beklemek nekadar doğru?

  • FxDev

    Cevap biraz geç oldu, kusura bakmayın.
    – İlk kısımdaki divider tanımlaması kurulacak timer için gerekli değerdir. Bunu datasheetten ilgili timer birimine bakarak anlayabilirsiniz.
    – İkincisi ise offset ile alakalı. Bildiğiniz gibi sinüs sinyalinin artı ve eksi alternansları vardır. Fakat PIC eksi alternansı ölçemeyeceği için bir offset değeri verilir. 3.3V/2=1.65V’luk bir offset ve max. 1.65V genlikli bir sinüs bileşeni PIC’in ADC’sine girdiğinde 1.65V sanki 0V’muş gibi davrandırılır. Bu gerilim de PIC’in içerisinde 512′ye denk gelmektedir. Dolayısıyla okunan değerden 512 çıkartılarak, asıl okunmak istenen değer – ve + şekilde programa işletilir.
    -ADC okuması kesmeden oldukça hızlı bir şekilde tamamlandığından herhangi bir sorun göremiyorum. Elbette normal uygulamalarda örneğin ADC_oku şeklinde bir flag set edilerek asıl ADC okuması main fonksiyonunda yapılmalıydı.

  • Overlord

    Merhaba FxDev öncelikle çalışma için tebrik ederim. Sinyalin frekansına göre LCD’deki şekillerin nasıl oluştuğunu anlamadım. Yani sütunlar nasıl o şekilde ekolayzır gibi davranıyo. Bunu anlatırsanız çok mutlu olurum.

  • FxDev

    Bunu birinci bölümde anlattım aslında. Uygulanan algoritma sayesinde sinyalin ufak bir bölümünden frekans bileşenleri ve bunların genlikleri ortaya çıkıyor. Bu da LCD’ye yansıtılıyor. Matematiksel işlemler ise biraz sayısal sinyal işleme bilgisi gerektiriyor.

  • Mustafa

    Frekans bandını yüksek frekanslara çektiğinizde neler oluyor. Aslında girişe müzik uygulasanız 10 veya 20 band spectrum olabilirmi Ayrıca LCD ile birlikte Led ilave edebilirmiyiz. Yani audio spectrum analyzer yapılabilirmi. Ayrıca tebrik ederim çok iyi bir çalışma .

  • FxDev

    Evet yapılabilinir örnek süresi artılırsa bence Goertzel bu işin altından FFT’ye göre rahatlıkla kalkabilir.

  • metin

    Etkili çalışmanız için sizi tebrik ederim Fırat bey, benim size bir sorum olacak cevaplaya bilirseniz sevinirim… Goertzel algoritması ile caller ID sinyali çözülebilir mi? bu örnekte siz sabit bir dalga şeklinin ( sabit derken sürekliliği kast ediyorum ) içerdiği frekansları yakalıyorsunuz. Ancak telefondan gelen Caller ID sinyalinde tarama işlemini yaparken (1200 ve 2200 hz) sinyalin diğer kısmını kaçırma gibi bir durum olabileceğini düşünüyorum. Yani caller ID için Goertzel kullanılabilir mi kullanılamaz mı paylaşırsanız sevinirim. İyi çalışmalar.

  • FxDev

    @Metin: Daha önce Caller ID sinyalleri ile uğraşmamış biri olarak cevabım teoriksel olacaktır. Yüksek bir örnekleme hızı ve fazla örnek sayısı alma ile bence Caller ID sinyalleri çözülebilir. Sonuçta algoritma size hangi sinyallerin baskın olduğunu söylüyor. Sonuç olarak hangi sayının geldiğini çözmek biraz da programa kalmış. Yalnız Caller ID başta data şeklinde geliyorsa yani binary 1-0 şeklinde bunu, bu algoritma ile çözmek imkansız. Bu algoritma bahsettiğim gibi baskın frekans çözümünde etkili. Caller ID data ise çözüm olanaksız.

  • metin

    Cevabınız için teşekkür ederim…

  • Ferhat YOL

    Merhabalar Fırat Bey.
    Paylaşımınız için Teşekkür ederim.
    Anlamadığım bir nokta var. İzninizle size danışmak istiyorum.
    Goertzel algooritmasında ADC örnekleme frekansı belirlenmesi istenen frekansın en az iki katı olması gerektiğini belirtmiştiniz.

    ADC örnekleme Frekansını nasıl belirlediniz. Daha doğru ADC ile örnek alma işlemini nasıl yaptınız. Burada Timer’in yaptığı işlem nedir?

    Birde iyi bir sonuç alabilmek için min. kaç adet örnek alınmalı?

    Bunları Açıklayabilirmisiniz.

  • FxDev

    @Ferhat: Örneğin siz 50-150-250-350Hz frekansları görmek istiyorsanız örnekleme frekansınız mutlaka 700Hz’den fazla olmalıdır. Bu bir kriter.

    Timer ise benim almak istediğim örnekleme frekansını oluşturuyor burada. Örnek verdiğim proje için örnekleme frekansında yani 700Hz’e set edilmiş bir Timer’ım var ve her Timer kesmesinde ADC kanalından bir örnek alıyorum.
    700Hz ile 10 kez örnek alırsanız çözünürlüğünüz 70Hz olur ve 50-150-250-350 bandındaki frekansları katları olmadığı için yakalayamazsınız. Bunun için 140 örnek aldığınızda ise 700:140=5Hz çözünürlüğünüz olacaktır ki 50 ve katlarını bu şekilde tanımlayabilirsiniz.

  • Ferhat YOL

    Cevabınız için Teşekkür Ederim.
    Biliyorum Bu iş için Pek uygun değiller ama Ben 18F serisi Bir Pic İle bu işlemi protonda yapmayı planlıyorum. Protonun Sample klasöründe FFT konusunda bir örnekde var.

    Benim amacım ilk başta en azından ses sinyalini 8 banda bölmek istiyorum.

    Bu iş için örnekleme frekansımı yani timer kesmesi frekansımı 1Khz olarak düşündüm.

    Kesme ile 100 defa örnek aldıktan sonra kesmeyi kapatıp aldığım örnekleri goertzel işlemlerine tabi tutacağım. Goertzel işlemleri bittip istediğim frekansı elde ettiğimde tekrar kesmeyi açıp örnek alıp bu sefer diğer görmek istediğim frekans ile goertzel işlemlerini başlatacağım. Aklımda böyle bir döngü var.

    Bu şekilde yapılan işlem sırası doğrumudur?

  • FxDev

    @Ferhat: 18F542 ile 128 noktalı FFT algoritması işletmiştim grafik LCD ile birlikte. Dikkat et FFT diyorum Goertzel değil. Bu örneği sitemde bulabilirsin.

    Dediğin sırada herhangi bir problem gözükmüyor.

  • Ferhat YOL

    Anladım Çok Teşekkür Ederim
    İyi Günler

  • Hakan K.

    Merhaba Fırat tebrik ederim güzel bir çalışma. Acaba HI-TECH demi yazdın kodları.

  • FxDev

    @Hakan: 24F serisinden bir PIC kullanıldığından kodlar C30 ile yazılmıştır.

  • Ali Deniz

    paylaşımın için tebrik ederim güzel paylaşım.18 f ler için olanından yapabilme şansınız var mı?c kodu lazım.Şebeke frekansı ölçse yeterli olacak.Veya şebeke frekansını analog kanaldan okumak için bir bilgin yöntemin var mı?