Generátor harmonického signálu

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

Zadanie

Na generovanie harmonického signálu nepoužite funkcie sin(), resp. cos(), ale vytvorte oscilátor (v reálnom čase) ako prenosovú funkciu 1/((s*T)^2 + 1). Úlohou je zmerať jeden bod frekvenčnej charakteristiky systému 1/(s*T_osc+1) na frekvencii omega = 1/T. T = 0,5 sek. Výstupom je signál s parametrami: A_0 + A_1*sin(omega*t + fi), kde A_0 = 128 a A_1 = 100.

Úlohou bolo generovať harmonický signal bez použitia funkcií sin() alebo cos(). Na tento účel bol použitý oscilátor realizovaný ako prenosová funkcia:

H(s) = 1 / ((s · T)^2 + 1)

Zároveň bolo potrebné zmerať jeden bod frekvenčnej charakteristiky systému:

H(s) = 1 / (s · T + 1)

pre frekvenciu: ω = 1 / T, kde T = 0.5 s.

Výstupný signal má mať tvar:

A₀ + A₁ · sin(ωt + φ)

kde: A₀ = 128, A₁ = 100.

Literatúra:

• Dokumentácia k doske Acrob
• Harmonic oscillator
• Transfer function
• Euler method

Analýza a opis riešenia

Cieľom riešenia je vytvoriť sinusový signal bez použitia matematických funkcií sin() alebo cos(). Tento problém sa rieši pomocou diskretneho oscilátora, ktorý vychádza z diferenciálnej rovnice harmonického kmitania.

Základom je rovnica:

y + ω²y = 0

Táto rovnica popisuje harmonické kmity a jej riešením sú funkcie sin() a cos(). To znamená, že ak vieme túto rovnicu numericky riešiť, vieme generovať sinus.

Teoretický základ a odvodenia

Laplaceova transformacia

Vychadzajme z prenosovej funkcie systému:

H(s) = Y(s) / X(s) = 1 / (T · s + 1)

Po úprave:

Y(s) · (T · s + 1) = X(s)

T · s · Y(s) + Y(s) = X(s)

Použitím inverznej Laplaceovej transformacie (kde s predstavuje deriváciu) dostaneme:

T · dy(t)/dt + y(t) = x(t)

Týmto získame diferenciálnu rovnicu systému v časovej oblasti.

Diskretizacia

Mikrokontrolér pracuje v diskretnom čase, preto je potrebné nahradiť derivácie rozdielmi medzi vzorkami.

Pre druhú deriváciu použijeme aproximáciu:

y ≈ (y[n] − 2y[n−1] + y[n−2]) / T_s²

Po dosadení do y + ω²y = 0:

y[n] = (2 / (1 + ω²T_s²)) · y[n−1] − (1 / (1 + ω²T_s²)) · y[n−2]

Táto rovnica predstavuje numerickú aproximáciu oscilátora, ale nie je ideálna z hľadiska stability amplitúdy.

Diskretny oscilátor (presný model)

Presnejší prístup vychádza z trigonometrických identít:

sin(A + B) = sin(A)cos(B) + cos(A)sin(B)

sin(A − B) = sin(A)cos(B) − cos(A)sin(B)

Po sčítaní:

sin(A + B) + sin(A − B) = 2 · sin(A) · cos(B)

Dosadením:

A = (n−1)θ, B = θ

dostaneme:

sin(nθ) + sin((n−2)θ) = 2 · cos(θ) · sin((n−1)θ)

Označením:

y[n] = sin(nθ)

vznikne rekurentný vzťah:

y[n] = 2 · cos(θ) · y[n−1] − y[n−2]

Tento vzťah generuje stabilný sinusový signal bez zmeny amplitúdy.

Výpočet parametrov

Platí:

θ = ω · T_s

kde: T = 0.5 s ω = 1 / T = 2 rad/s SAMPLE_RATE = 1000 Hz T_s = 0.001 s

θ = 0.002

Aproximacia cos()

Keďže nie je dovolené použiť funkciu cos(), použije sa Taylorov rozvoj:

cos(θ) ≈ 1 − θ²/2

Z toho:

2cos(θ) ≈ 2 · (1 − θ²/2)

Inicializácia oscilátora

Pre správnu činnosť oscilátora sú potrebné počiatočné hodnoty:

y1 = 1, y2 = 1 − θ²/2

Generovanie signalu

Oscilátor generuje hodnoty v rozsahu:

[-1, 1]

Požadovaný výstup:

x = A0 + A1 · y

kde:

A0 = 128 A1 = 100

Diskretizacia systému 1 / (sT + 1)

Zo spojitej rovnice:

T · dy/dt + y = x

Použitím Eulerovej metódy:

dy/dt ≈ (y[n] − y[n−1]) / T_s

Po úprave:

y[n] = x − (T / T_s) · y[n] − (T / T_s) · y[n−1]

y[n] · (1 + (T / T_s)) = x − (T / T_s) · y[n−1]

y[n] = (T_s / (T + T_s)) · x + (T / (T + T_s)) · y[n−1]

Po zavedení:

α = T_s / (T + T_s)

dostaneme praktický tvar:

y[n] = α · x + (1 − α) · y[n−1]

y[n] = y[n−1] + α · (x − y[n−1])

Algoritmus a program

Algoritmus programu využíva diskretny oscilátor a numericku aproximáciu systému 1 / (sT + 1). Základné výpočty prebiehajú v prerušení Timer1 s frekvenciou 1 kHz:

• výpočet oscilátora
• generovanie vstupu x
• výpočet výstupu systému
• výstup cez PWM
• odosielanie dát cez UART

main.cuart.huart.c

#define F_CPU 16000000 UL

#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) {}

}

#define set_bit(ADDRESS, BIT)(ADDRESS |= (1 << BIT))
#define clear_bit(ADDRESS, BIT)(ADDRESS &= ~(1 << BIT))

#ifndef UART_H_
#define UART_H_

#include <stdio.h>

#define BAUD_PRESCALE(((F_CPU / (BAUDRATE * 16 UL))) - 1)

void uart_init(void);
int uart_putc(char c, FILE * stream);
void uart_puts(const char * s);
char uart_getc(void);
void delay(int delay);

#endif /* UART_H_ */

#include <avr/io.h>

#include <util/delay.h>

#include "uart.h"

void uart_init(void) {
    #include <util/setbaud.h>

    UBRR0H = UBRRH_VALUE;
    UBRR0L = UBRRL_VALUE;
    #if USE_2X
    UCSR0A |= (1 << U2X0);
    #else
    UCSR0A &= ~(1 << U2X0);
    #endif

    UCSR0C = _BV(UCSZ01) | _BV(UCSZ00);
    UCSR0B = _BV(RXEN0) | _BV(TXEN0);
}

int uart_putc(char c, FILE * stream) {
    if (c == '\n') {
        loop_until_bit_is_set(UCSR0A, UDRE0);
        UDR0 = '\r';
    }

    loop_until_bit_is_set(UCSR0A, UDRE0);
    UDR0 = c;

    return 0;
}

void uart_puts(const char * s) {}

char uart_getc(void) {
    loop_until_bit_is_set(UCSR0A, RXC0);
    return UDR0;
}

void delay(int delay) {
    for (int i = 1; i <= delay; i++)
        _delay_ms(1);
}

Zdrojový kód: MykytaOleksandr_sources.zip

Overenie

Funkcia systému bola overená pomocou výpisu dát cez UART. Do sériového portu sa posielajú dvojice hodnôt:

x, y_sys

Tieto hodnoty je možné zobraziť napríklad pomocou SerialPlot, kde je viditeľný vstupný sinusový signal a výstup systému.

Namerené priebehy vstupu a výstupu zo Serial Plot

• 
  x_max = 227
• 
  x_min = 28
• 
  y_max = 198
• 
  y_min = 57

Pri overovaní bolo sledované, či:

• amplitúda výstupu zodpovedá očakávaniu
• systém vykazuje fázový posun
• priebeh signalu je stabilný bez driftu
• výstup reaguje správne na vstupný sinus

Zo signalov je možné pozorovať zmenu amplitúdy a fázový posun, čo predstavuje bod frekvenčnej charakteristiky systému.

Video:

Video demonštruje reálne správanie systému. Je na ňom viditeľné generovanie sinusového signalu, jeho spracovanie systémom a výstupné dáta zobrazené v SerialPlot. Video slúži na overenie, že implementácia funguje podľa očakávania.

Čo by som urobil inak

Na ďalšíkrát by som si lepšie rozmyslel, ako presne riešiť diskretizáciu, lebo ten jednoduchý prístup funguje, ale nie je úplne ideálny.

A ešte by som si skôr urobil priamo meranie aj na mikrokontroléri v kode, aby som nemusel všetko riešiť cez SerialPlot.