Widma częstotliwościowe prostych sygnałów

Tworzenie aplikacji analizatora widma najlepiej rozpocząć od wykonania czegoś prostego. Ponieważ interesuje nas widmo częstotliwościowe to warto sprawdzić jak wygląda ono dla kilku prostych sygnałów. W tym celu możemy wygenerować cztery sygnały, które zawierają od jednej do czterech składowych harmonicznych:

if(mode == 1){
  for(i=0; i<512; ++i){
    buffer_input[i] = (float32_t) 30*sin(2*PI*freq*i*dt);
  }
}
else if(mode == 2){
  for(i=0; i<512; ++i){
    buffer_input[i] = (float32_t) 20*sin(2*PI*freq*i*dt) +
                                  10*sin(2*PI*2*freq*i*dt);
  }
}
else if(mode == 3){
  for(i=0; i<512; ++i){
    buffer_input[i] = (float32_t) 15*sin(2*PI*freq*i*dt) +
                                  10*sin(2*PI*2*freq*i*dt) +
                                  5*sin(2*PI*3*freq*i*dt);
  }
}
else if(mode == 4){
  for(i=0; i<512; ++i){
    buffer_input[i] = (float32_t) 15*sin(2*PI*freq*i*dt) +
                                  10*sin(2*PI*2*freq*i*dt) +
                                  3*sin(2*PI*3*freq*i*dt) +
                                  2*sin(2*PI*4*freq*i*dt);
  }
}

Generowanych jest 512 próbek sygnału i liczba ta wynika z wymagań narzuconych przez funkcję wyznaczającą widmo i chęci przedstawienia jak najdłuższego przebiegu na wyświetlaczu (wykorzystane będą tylko pierwsze 256). Wartość zmiennej freq w toku działania aplikacji demonstracyjnej jest zmieniana od 1 do 1024, a zmienna dt jest wirtualnym okresem próbkowania, który wynosi 0.001. Dalej można wyświetlić wygenerowane sygnały:

// Usunięcie poprzedniego przebiegu sygnału
LCD_SetForegroundColor(Black);
for(i=0; i<255; ++i){
  x_LCD_DrawLine_2(i + 32, (uint16_t)(buffer_input_copy[i] + 50),
                   i + 33, (uint16_t)(buffer_input_copy[i+1] + 50));
}
// Rysowanie nowego przebiegu sygnału
LCD_SetForegroundColor(Yellow);
for(i=0; i<255; ++i){
  x_LCD_DrawLine_2(i + 32, (uint16_t)(buffer_input[i] + 50),
                   i + 33, (uint16_t)(buffer_input[i+1] + 50));
}

// Wyświetlenie informacji o częstotliwościach składowych harmonicznych sygnału
if(mode == 1){
  sprintf(text,"F=%4d [1]",freq);
  x_LCD_DrawText(32, 100, text, LeftJustify, NormalOrientation);
}
else if(mode == 2){
  sprintf(text,"F1=%4d [1], F2=%4d [1]",freq,freq*2);
  x_LCD_DrawText(32, 100, text, LeftJustify, NormalOrientation);
}
else if(mode == 3){
  sprintf(text,"F1=%4d [1], F2=%4d [1]",freq,freq*2);
  x_LCD_DrawText(32, 100, text, LeftJustify, NormalOrientation);
  sprintf(text,"F3=%4d [1]",freq*3);
  x_LCD_DrawText(32, 110, text, LeftJustify, NormalOrientation);
}
else if(mode == 4){
  sprintf(text,"F1=%4d [1], F2=%4d [1]",freq,freq*2);
  x_LCD_DrawText(32, 100, text, LeftJustify, NormalOrientation);
  sprintf(text,"F3=%4d [1], F4=%4d [1]",freq*3,freq*4);
  x_LCD_DrawText(32, 110, text, LeftJustify, NormalOrientation);
}

// Zachowanie kopii pierwszej połowy próbek aktualnego sygnału
for(i=0; i<256; ++i){
	buffer_input_copy[i] = buffer_input[i];
	buffer_output_mag_copy[i] = buffer_output_mag[i];
}

Mając przygotowany sygnał wykorzystujemy trzy funkcje z biblioteki DSP i wyznaczamy najpierw transformatę Fouriera danego sygnału, a następnie obliczamy wartości modułów dla każdej pary liczb rzeczywistej i urojonej. Na koniec możemy wyszukać składową sygnału (częstotliwość), która ma największą amplitudę i przeskalować tablicę modułów do zakresu nadającego się do wyświetlenia:

// Wyznaczenie transformaty Fouriera
arm_rfft_f32(&S, buffer_input, buffer_output);
// Obliczenie modułów
arm_cmplx_mag_f32(buffer_output, buffer_output_mag, 512);
// Znalezienie składowej harmonicznej sygnału o największej amplitudzie
arm_max_f32(buffer_output_mag, 512, &maxvalue, &maxvalueindex);
// Skalowanie wartości modułów
for(i=0; i<512; ++i){
	buffer_output_mag[i] = 100*buffer_output_mag[i]/maxvalue;
}

Pierwsza funkcja przyjmuje trzy argumenty: wskaźnik do struktury typu arm_rfft_instance_f32 opisującej sposób działania transformaty FFT, wskaźnik do buforu wejściowego z wartościami próbek sygnału oraz wskaźnik do buforu wyjściowego w którym będą zapisane zespolone wartości transformaty sygnału. Ogólną ideę funkcjonowania tej funkcji pokazano na rysunku 9. Ważne jest aby w buforze wejściowym były zapisane rzeczywiste wartości (UWAGA: zostaną one zmodyfikowane po wykonaniu się funkcji), czyli:

{real[0], real[1], real[2], real[3], ..}

Natomiast w buforze wyjściowym, który powinien być dwa razy większy od wejściowego, oprócz wartości rzeczywistych będą znajdować się również wartości urojone:

{real[0], imag[0], Real[1], imag[1], ...}

W przypadku struktury przekazywanej w pierwszym argumencie nie ma potrzeby ręcznego uzupełniania jej, gdyż można do tego celu wykorzystać funkcję arm_rfft_init_f32():

// Wystarczą deklaracje dwóch struktur
arm_rfft_instance_f32 S;
arm_cfft_radix4_instance_f32 S_CFFT;

// oraz inicjalizacja, gdzie:
//  arg#3 – liczba próbek w buforze wejściowym; możliwe wartości to: 128, 512, 2048
//  arg#4 – kierunek transformacji: 0 – prosta, 1 – odwrotna
//  arg#5 – uporządkowanie wartości w buforze wyjściowym: 0 – odwrócone, 1 - normalne
arm_rfft_init_f32(&S, &S_CFFT, 512, 0, 1);

Rys. 9. Idea działania funkcji arm_rfft_f32()

Druga funkcja przyjmuje również trzy argumenty: wskaźnik do buforu wejściowego z wartościami par liczb rzeczywistych i urojonych, wskaźnik do buforu wyjściowego do którego zostaną zapisane wartości rzeczywiste oraz ostatni parametr – ilość liczb zespolonych. Funkcja oblicza moduł dla każdej pary wartości liczby rzeczywistej i urojonej zgodnie ze wzorem:

Trzecia funkcja znajduje maksymalną wartość w tablicy buffer_output_mag o długości 512 elementów i zapisuje ją do zmiennej maxvalue, a indeks komórki w jakiej się ona znajduje jest zapisywany w zmiennej maxvalueindex.

Fot. 10. Sygnał o jednej składowej harmonicznej i jego widmo

Fot. 11. Sygnał o dwóch składowych harmonicznych i jego widmo

Fot. 12. Sygnał o trzech składowych harmonicznych i jego widmo

Fot. 13. Sygnał o czterech składowych harmonicznych i jego widmo

Na koniec, po przeskalowaniu wartości modułów, można wyświetlić otrzymane wyniki (fotografie 10…13). Wystarczy przedstawić pierwsze 256 elementów tablicy buffer_output_mag, gdyż druga połowa stanowi lustrzane odbicie pierwszej:

// Usunięcie poprzedniego widma
LCD_SetForegroundColor(Black);
for(i=0; i<256; ++i){
  x_LCD_DrawLine_2(i + 32,
                   130 + 100,
                   i + 32,
                   130 + (uint16_t)(100 - (uint16_t)buffer_output_mag_copy[i+1]));
}
// Rysowanie nowego widma
LCD_SetForegroundColor(Green);
for(i=0; i<256; ++i){
  x_LCD_DrawLine_2(i + 32,
                   130 + 100,
                   i + 32,
                   130 + (uint16_t)(100 - (uint16_t)buffer_output_mag[i+1]));
}