Generátor harmonického signálu

Záverečný projekt predmetu MIPS / LS2026 - Oleksandr Mykyta

Zadanie

Ulohou je generovat harmonicky signal bez pouzitia funkcii sin() alebo cos(). Namiesto toho sa ma pouzit oscilator realizovany ako prenosova funkcia:

$$H(s) = \frac{1}{(s \cdot T)^2 + 1}$$

Zaroven je potrebne zmerat jeden bod frekvencnej charakteristiky systemu:

$$H(s) = \frac{1}{s \cdot T + 1}$$

na frekvencii ω = 1 / T, pricom T = 0,5 s.

Vystupny signal ma mat tvar:

$$A_0 + A_1 \cdot \sin(\omega t + \varphi)$$

kde: A₀ = 128 A₁ = 100

Analýza a opis riešenia

Cielom riesenia je vytvorit sinusovy signal bez pouzitia matematickych funkcii sin() alebo cos(). Tento problem sa riesi pomocou diskretneho oscilatora, ktory vychadza z diferencialnej rovnice harmonickeho kmitania.

1. 1. 1. Laplaceova transformacia a prenosova funkcia

Vychadzajme z nasej prenosovej funkcie H(s), co je v podstate len pomer vystupu k vstupu v Laplaceovej oblasti:

$$H(s) = \frac{Y(s)}{X(s)} = \frac{1}{T \cdot s + 1}$$

Ak chceme vediet, ako sa to sprava v case, potrebujeme to transformovat naspat. Vynasobime rovnicu krizom:

$$Y(s) \cdot (T \cdot s + 1) = X(s)$$ $$T \cdot s \cdot Y(s) + Y(s) = X(s)$$

Teraz pouzijeme inverznu Laplaceovu transformaciu. Hlavne pravidlo tu je, ze premenna s sa v casovej domene sprava ako derivacia: $\mathcal{L}^{-1}\{s \cdot Y(s)\} = \frac{dy}{dt}$. Po aplikovani inverznej transformacie dostaneme diferencialnu rovnicu:

$$T \cdot \frac{dy(t)}{dt} + y(t) = x(t)$$

Tymto sme si potvrdili, ze nasa prenosova funkcia v Laplaceovi a diferencialna rovnica v case su vlastne len dva pohlady na to iste.

1. 1. 1. Diskretizácia

Mikrokontroler pracuje v diskretnom case, preto je potrebne nahradit derivacie rozdielmi medzi vzorkami (metoda konecnych diferencii).

Pouzie sa aproximacia druhej derivacie pre oscilator:

$$y \approx \frac{y[n] - 2y[n-1] + y[n-2]}{T_s^2}$$

1. 1. 1. Diskrétny oscilátor

Na generovanie stabilneho sinusu sa pouzije presny diskretny model zalozeny na trigonometrickej identite:

$$\sin(n\theta) = 2\cos(\theta)\cdot\sin((n-1)\theta) - \sin((n-2)\theta)$$

Po oznaceni:

$$y[n] = \sin(n\theta)$$

dostaneme rekurentny vztah:

$$y[n] = 2\cos(\theta)\cdot y[n-1] - y[n-2]$$

Tento vztah generuje stabilny sinusovy signal bez zmeny amplitudy. Pre vypocet parametrov plati $\theta = \omega \cdot T_s$.

1. 1. 1. Aproximácia cos()

Kezde nie je dovolene pouzit funkciu cos(), pouzie sa aproximacia (Taylorov rozvoj):

$$\cos(\theta) \approx 1 - \frac{\theta^2}{2}$$

Tento vyraz sa pouzie ako koeficient oscilatora v programe.

1. 1. 1. Systém 1 / (sT + 1) a koeficient alpha

Aby mikrokontroler vedel vyuzit system 1 / (sT + 1), musime ho previest do diskretneho casu pomocou Eulerovej metody. Derivaciu nahradime diferenciami:

$$\frac{dy}{dt} \approx \frac{y[n] - y[n-1]}{T_s}$$

Ked to dosadime do diferencialnej rovnice:

$$T \cdot \frac{y[n] - y[n-1]}{T_s} + y[n] = x$$

Po uprave dostaneme:

$$y[n] = \left(\frac{T_s}{T + T_s}\right) \cdot x + \left(\frac{T}{T + T_s}\right) \cdot y[n-1]$$

Ak si oznacime α = Ts / (T + Ts), tak to vyzera ovela krajsie:

$$y[n] = \alpha \cdot x + (1 - \alpha) \cdot y[n-1]$$

Co je to iste ako:

$$y[n] = y[n-1] + \alpha \cdot (x - y[n-1])$$

Pre dane hodnoty α ≈ 0,001996.

Algoritmus a program

Algoritmus programu vyuziva Timer1 prerusenie, v ktorom sa vypocitava oscilator a nasledne filter. Vystup sa posiela do PWM modulu a data cez UART.

AVR C-code

#define F_CPU 16000000UL

#include <avr/io.h>
#include <avr/interrupt.h>
#include <stdio.h>
#include "uart.h"

#define SAMPLE_RATE 1000.0
#define T 0.5

#define A0 128
#define A1 100

float OSC_COEFF;

volatile float y = 0;
volatile float y1 = 0;
volatile float y2 = 0;

volatile float y_sys = 0;

float alpha;

FILE mystdout = FDEV_SETUP_STREAM(uart_putc, NULL, _FDEV_SETUP_WRITE);

ISR(TIMER1_COMPA_vect)
{
    y = OSC_COEFF * y1 - y2;

    y2 = y1;
    y1 = y;

    float x = A0 + A1 * y;

    y_sys = y_sys + alpha * (x - y_sys);

    OCR0A = (uint8_t)(y_sys);

    printf("%d,%d\n", (int)x, (int)y_sys);
}

void timer1_init()
{
    TCCR1B |= (1 << WGM12);

    OCR1A = 15999;

    TCCR1B |= (1 << CS10);

    TIMSK1 |= (1 << OCIE1A);
}

void pwm_init()
{
    DDRD |= (1 << PD6);

    TCCR0A |= (1 << COM0A1) | (1 << WGM01) | (1 << WGM00);
    TCCR0B |= (1 << CS01);
}

int main(void)
{
    uart_init();
    stdout = &mystdout;

    pwm_init();
    timer1_init();

    float Ts = 1.0 / SAMPLE_RATE;

    alpha = Ts / (T + Ts);

    float theta = (1.0 / T) * (1.0 / SAMPLE_RATE);

    y1 = 1.0;
    y2 = 1.0 - (theta * theta) / 2.0;

    OSC_COEFF = 2.0 * (1.0 - (theta * theta) / 2.0);

    sei();

    while (1)
    {
    }
}

Overenie

Funkcia systemu bola overena pomocou vypisu dat cez UART. Do serioveho portu sa posielaju dvojice hodnot: x, y_sys. Tieto hodnoty je mozne zobrazit napriklad pomocou Serial Plotteru, kde je viditelny vstupny sinusovy signal a vystup systemu.

Zo signalov je mozne pozorovat zmenu amplitudy a fazovy posun, co predstavuje bod frekvencnej charakteristiky systemu. System sa chova presne tak, ako hovori teoria.

Čo by som urobil inak

Pri rieseni by bolo mozne pouzit presnejsiu metodu diskretizacie systemu 1/(sT + 1), napriklad bilinearnu transformaciu, ktora by zlepsila presnost modelu.

Taktiez by bolo mozne implementovat presnejsi vypocet cos(θ) bez aproximacie, napriklad pomocou lookup tabulky.

Dalsim zlepsenim by mohlo byt automaticke vyhodnotenie amplitudy a fazoveho posunu priamo v mikrokontroleri namiesto spracovania na PC.