Тема: Дискретне перетворення Фур’є. Розрахунок коефіцієнтів ДПФ та реалізація на ЦСП С55х.

 

Дискретне перетворення Фур’є (ДПФ) є однією з найпоширеніших функцій цифрової обробки сигналів. Для його реалізації у пакеті MATLAB є спеціальна функція – dft.m. Алгоритм ДПФ можна реалізувати за допомогою програм, написаних на інших мовах для програмування. Наприклад, нижче наведено текст такої функції, написаний мовою С, що реалізує виконання алгоритму ДПФ для 128 точок:

 

/* Experiment A – expa.c */

#include “input.dat”

#define N 128

extern void dft_128(int *, int *);

extern void mag_128(int *, int *);

int Xin[2*N];

int Xout[2*N];

int Spectrum[N];

void main()

{

int i,j;

for (j=0,i=0;i<N;i++)

{

Xin[j++] = input[i]; /* Get real sample */

Xin[j++] = 0; /* Imaginary sample=0 */

}

dft_128(Xin, Xout);     /* DFT routine */

mag_128(Xout, Spectrum); /* Compute spectrum */

}

 

Обчислення ДПФ вимагає організації вкладених циклів, виконання операції перемноження вибірок комплексних чисел та створення масивів комплексних коефіцієнтів перетворення. У даній лабораторній роботі необхідно написати підпрограму мовою асемблер для реалізації алгоритму розрахунку ДПФ.

    За умовою, вибірка вхідного сигналу представляє собою комплексне число: . Коефіцієнт перетворення Фур’є також є комплексним і виглядає наступним чином: . Добуток двох комплексних чисел x(n) та є комплексне число, яке представляється так:

,            (1)

де нижні індекси i та r вказують на уявну та дійсну частини комплексних чисел. Рівняння (1) можна записати у іншому вигляді:

                    (2)

для n=0, 1… N-1 де:

                    (3а)

                    (3б)

    У програмі expA.c, що наведена вище, використовуються два масиви Xin[2*N] та Xout[2*N] для зберігання вхідних та вихідних даних у вигляді комплексних чисел (дійсної та уявної частин). Вхідні дані для завдання генеруються за допомогою програми пакету MATLAB, текст якої наведено нижче:

 

fs =8000;                    %Sampling frequency in Hz

f1=500;                    %1st sinewave frequency in Hz

f2=1000;                    %2st sinewave frequency in Hz

f3=2000;                    %3st sinewave frequency in Hz

n=[0:127];                    %n=0, 1, …, 127

w1=2*pi*f1/fs;

w2=2*pi*f2/fs;

w3=2*pi*f3/fs;

x=0.3*sin(w1*n)+0.3*sin(w2*n)+ 0.3*sin(w3*n);

fid=fopen(‘input.dat’,’w’);        %Open input.dat for write

fprintf(fid, ‘int input[128]={\n’);

fprintf(fid, ‘%5d,\n’, round(x(1:127)*32767));

fprintf(fid, ‘%5d};\n’, round(x(128)*32767));

fclose(fid);                    %close file input.dat

 

Ця програма генерує 128 відліків вхідного сигналу, які зберігаються у файлі input.dat (програма знаходиться у директорії вхідних даних до завдання). Файл input.dat зчитується програмою expA.c. Перед початком роботи з середовищем CCS необхідно відкрити файл in4.m у пакеті MATLAB, та створити за його допомогою файл input.dat у робочій директорії проекту.

 

Завдання А

Розрахунок коефіцієнтів перетворення Фур’є

 

Для розрахунку коефіцієнтів перетворення Фур’є необхідно мати набір вибірок синусної та косинусної функцій. Його можливо зробити у вигляді табличної функції, або за допомогою спеціальної програми. В даному завданні буде використовуватися спеціальна програма – синусно-косинусний генератор, представлений файлом SIN_COS.asm. До функції sine_cos(angle, Wn) передаються тільки два аргументи. Переший аргумент angle – вхідний кут у радіанах, передається за допомогою регістру зберігання тимчасового значення Т0. Другий аргумент передається через допоміжний регістр AR0 і використовується як вказівник області пам’яті даних,
де зберігаються результати розрахунку Wn ( дві комірки пам’яті – для синусної sin(angle) та косинусної cos(angle) складових). Цю підпрограму на мові асемблер можна викликати як функцію з програми на мові С. Наведений нижче приклад демонструє виклик підпрограми генератора sine_cos програмою на мові С (вкладені цикли використовуються для генерації коефіцієнтів перетворення Фур’є):

 

 

#define N 128

#define TWOPIN 0×7FFF >>6        /*2pkn/N, N=128*/

int n, k, angle;

int Wn[2];                /*Wn[0]=cos(angle), Wn[1]=sin(angle)*/

for(k=0; k<N; k++)

{

    for(n=0; n<N; n++)

    {

        angle=TWOPIN*k*n;

        sine_cos(angle, Wn);

    }

}

Приклад фрагменту програми на мові асемблер, яка викликає підпрограму sine_cos, наведений нижче:

 

mov    *(angle), T0            ;Pass the first argument in T0

amov    #Wn,XAR0        ;Pass the second argument in AR0

call sine_cos            ;Call sine_cos subroutine

 

    Якщо програма розрахунку коефіцієнтів перетворення Фур’є написана мовою асемблер, то для реалізації вкладених циклів використовуються оператори повтору. Вкладений цикл цієї програми містить декілька різних операторів, а для їх циклічного виконання використовується оператор повтору блоку команд (rptb). Він застосовується для реалізації як вкладеного, так і зовнішнього циклів. Процесор С55х має два лічильники циклів – регістри BRC0 та BRC1. Коли реалізується вкладений цикл, як лічильник вкладеного циклу повинен використовуватись регістр BRC1, а для зовнішнього циклу використовується регістр BRC0. Процесор С55х дозволяє автоматично перезавантажувати лічильник вкладеного циклу BRC1, коли значення лічильника зовнішнього циклу змінюється. Нижче наведено приклад використання регістрів BRC0 та BRC1 у вкладених циклах (кількість повтору блоку команд – N разів):

 

mov #N–1, BRC0

mov #N–1, BRC1

rptb outer_loop–1

(more outer loop instructions…)

rptb inner_loop–1

(more inner loop instructions…)

inner_loop

(more outer loop instructions…)

outer_loop

    Значення коефіцієнта перетворення Фур’є є функція від кута, який в свою чергу залежить від двох змінних n (номер гармоніки) та k (фазовий зсув) таким чином:

angle=(2p/N)kn.                        (4)

При використанні операцій з фіксованою крапкою, для розширення динамічного діапазону, значення p в програмі синусно-косинусного генератора подається шістнадцятирічним числом 0×7FFF. Тому, значення кута, яке використовується для обчислення коефіцієнтів перетворення Фур’є при N=128, може бути записано так:

angle=(2*0×7FFF/128)*k*n=(0×1FF) *k*n.            (5)

Коефіцієнти n=0, 1, …, N–1 формують вкладений цикл, а коефіцієнти k=0, 1, …, N–1 зовнішній цикл. Нижче наведено текст підпрограми, написаної мовою асемблер, яка обчислює спектральні складові ДПФ для N=128:

 

;

; DFT_128 – compute 128 points DFT

;

; Entry T0: AR0: pointer to complex input sample buffer

; AR1: pointer to complex output sample buffer

; Return: None

;

.def    _dft_128

.ref    _sine_cos    

 

N .set 128

TWOPIN .set 0x7fff>>6 ; 2*PI/N, N=128

.bss Wn,2 ; Wn[0]=Wr, Wn[1]=Wi

.bss angle,1 ; Angle for sine-cosine function

.text    

 

_dft_128

pshboth XAR5 ; Save AR5

bset SATD

mov #N-1,BRC0 ; Repeat counter for loop N times

mov #N-1,BRC1 ; Repeat counter for loop N times

mov XAR0,XAR5 ; AR5 pointer to sample buffer

mov XAR0,XAR3

mov #0,T2 ; k = T2 = 0

rptb outer_loop-1 ; for(k=0;k<N;k++) {

mov XAR3,XAR5 ; Reset x[] pointer

mov #TWOPIN<<#16,AC0; hi(AC0) = 2*PI/N

mpy T2,AC0

mov #0,AC2 ; Xr[k] = 0

mov #0,AC3 ; Xi[k] = 0

mov #0,*(angle)    

mov AC0,T3 ; angle=2*PI*k/N

rptb inner_loop-1 ; for(n=0;n<N;n++) {

mov *(angle),T0 ; T0=2*PI*k*n/N

mov *(angle),AC0

add T3,AC0

mov AC0,*(angle) ; Update angle

amov #Wn,XAR0 ; AR0 is the pointer to Wn

call _sine_cos ; sine_cos(angle, Wn)

bset SATD ; sine_cos turn off FRCT & SATD

macm40 *AR5+,*AR0,AC2 ; XR[k]+Xin[n]*Wr

macm40 *AR5-,*AR0+,AC3 ; XI[k]+Xin[n+1]*Wr

masm40 *AR5+,*AR0,AC3 ; XI[k]+Xin[n+1]*Wr-Xin[n]*Wi

macm40 *AR5+,*AR0-,AC2 ; XR[k]+Xin[n]*Wr+Xin[n+1]*Wi

inner_loop ; } end of inner loop

mov hi(AC2<<#-5),*AR1+

mov hi(AC3<<#-5),*AR1+

add #1,T2

outer_loop ; } end of outer loop

popboth XAR5

bclr SATD

ret

.end

 

Фрагмент тексту програми мовою асемблер, який призначений для обчислення рівняння (3), наведено нижче:

 

mov        #0, AC2

mov        #0, AC3

rptb        inner_loop-1

macm40    *AR5+, *AR0, AC2        ;Xin[n]*Wr

macm40    *AR5–, *AR0+, AC3        ; Xin[n+1]*Wr

macm40    *AR5+, *AR0, AC3        ;Xi[k]=Xin[n+1]*Wr–Xin[n]*Wi

macm40    *AR5+, *AR0–, AC2        ;Xr[k]=Xin[n]*Wr+Xin[n+1]*Wi

Inner_loop

Оскільки для обчислення однієї складової ДПФ необхідно додати N проміжних результатів, то повинно врахувати можливість виникнення переповнення у процесі виконання обчислень. Оператор macm40 дає можливість використовувати захисні біти акумулятора, що дозволяє виконувати операцію множення з накопиченням з розширенням результату до 40 біт.

 

Завдання Б

Реалізація ДПФ за допомогою процесора С55х

 

У цьому розділі буде розроблено та випробувано програму реалізації ДПФ для N=128. Головна програма expA.c на мові С, текст якої наведено вище, викликає підпрограму dft_128.asm, написану мовою асемблер, для обчислення складових частин ДПФ.

    Вхідний файл даних input.dat – це ASCII файл, що містить 128 дискретизованих відліків, отриманих з частотою дискретизації 8 КГц. Перш за все, головна програма заповнить нулями уявні частини в масиві комплексних вхідних даних Xin[2*N]. Потім викликається підпрограма dft_128.asm. Результат обчислення 128 вхідних вибірок комплексного сигналу зберігається у масиві вихідних значень Xout[2*N]. Потім, отримані значення надсилаються до підпрограми mag_128.asm, яка розраховує амплітудні значення спектральних складових вхідного сигналу. Отримані значення зберігаються у масиві Spectrum[N], який використовується для графічного відображення амплітудно-частотної характеристики. Текст підпрограми mag_128.asm на мові асемблер, наведено нижче:

 

;

; Compute the magnitude of the input complex sample

;

; Entry: AR0: input buffer pointer

; AR1: output buffer pointer

; Exit: None

;

.def _mag_128

N .set 128

 

.text

_mag_128

bset SATD

pshboth XAR5

mov #N-1,BRC0 ; Set BRC0 for loop N times

mov XAR0,XAR5

bset FRCT

rptblocal mag_loop-1

mpym *AR0+,*AR5+,AC0 ; Xr[i]*Xr[i]

macm *AR0+,*AR5+,AC0 ; Xr[i]*Xr[i]+Xi[i]*Xi[i]

mov hi(saturate(AC0)),*AR1+

mag_loop

popboth XAR5

bclr SATD

bclr FRCT

ret

.end

Для виконання завдання зробіть наступні кроки:

 

1. Створіть новий проект у середовищі CCS; назвіть його expА та збережіть його у відповідній директорії. Напишіть програму expА.с на основі наведеного вище коду та збережіть її у відповідній директорії. Скопіюйте командний файл лінкера exp4.cmd з директорії вхідних даних до завдання. Додайте до проекту файли expА.с, exp4.cmd, mag_128.asm, dft_128.asm та SIN_COS.asm. Підключить бібліотеку засобів динамічної підтримки rst55.lib (розташована у директорії C:\ti\c5500\cgtools\lib). Запустіть програму на компіляцію. Використайте згенерований файл input.dat (розміщений в робочій директорії), як вхідний сигнал для обробки.

2. Завантажте програму до процесора. Відкрийте графічне вікно для перегляду вхідного сигналу. Параметри налагодження графічного вікна наведені на рис.1.

Рис.1 Приклад настройки графічного вікна для виведення вхідного сигналу.

 

3. Відкрийте графічне вікно для перегляду спектру вхідного сигналу. Параметри налагодження графічного вікна наведені на рис.2.

 

Рис.2 Приклад настройки графічного вікна для виведення спектру вхідного сигналу.

 

4. Відкрийте графічне вікно для перегляду результатів ДПФ. Параметри налагодження графічного вікна наведені на рис.3.

 

Рис.3 Приклад настройки графічного вікна для виведення результатів ДПФ.

5. Запустіть програму на виконання за допомогою команди RUN. Поясніть отриманий результат.

 

 

Завдання для самостійної роботи.

 

А. Доопрацюйте вхідні дані в середовищі MATLAB таким чином, щоб в них були присутні інші частотні складові. Перевірте результат обробки даних за допомогою ДПФ.

Б. Доопрацюйте програму – виконайте обробку вхідного сигналу для іншої кількості відліків (кількість відліків повинна бути кратною ступеням числа 2).