Быстрое преобразование Фурье. Практика использования. Часть 2

Posted by bullvinkle under Журнал, Статьи

Традиционно, измерительный прибор представляет собой автономное специализированное устройство, к которому подключаются анализируемые входные сигналы и на выходе которого имеется некий результат. Причем внутренняя архитектура остается неизменной, в отличие от виртуальных приборов. И если для изменения функциональности в первом, требуется существенная переработка схемы и конструкции, то для вторых достаточно изменить программу. Продолжая наш цикл по практике использования быстрого преобразования Фурье [1], сегодня мы с вами научим наш виртуальный прибор работать не только с низкочастотной частью диапазона и аудиоустройствами, но и с высокоскоростной платой сбора данных и для удобства принимать команды управления основными режимами программы с пульта… Данная статья рассчитана в помощь программистам и инженерам-разработчикам в области цифровой обработки сигналов (DSP).

Быстрое преобразование Фурье. Практика использования. Часть 2

raxp www.programmersforum.ru

DSP – (Digital Signal Processing) преобразование сигналов, представленных в цифровой форме

АЦП – аналого–цифровой преобразователь

Зоны Найквиста – зеркальные отображения спектра при использовании частот выше половины частоты дискретизации (частоты Найквиста)

БПФ/ДПФ – быстрое / дискретное преобразование Фурье.

Зачем-же все это нужно? Как мы уже знаем, виртуальные приборы (ВП) благодаря гибкости в их построении все больше вытесняют дорогостоящие автономные аппаратные решения, таких как осциллографы, спектроанализаторы и др. При этом пользователь не ограничен в выборе средств для анализа и обработки информации, что сводится лишь к изменению программного обеспечения. Приведем пример: вы купили «имиджевый» осциллограф, используете его до поры до времени… и вот настает момент, когда меняется задача и вам понадобился анализ спектра АЦП. Снова затраты? Рассмотрим подробнее…

Краткий экскурс…

В настоящее время основополагающим принципом цифровой обработки сигналов (ЦОС) является преобразование аналогового сигнала в цифровой на промежуточной частоте (перенос спектра). Это позволяет исключить такие недостатки аналогового способа формирования квадратурных сигналов, как: невысокие стабильность и нелинейность, неидентичность каналов, смещение фазы и трудности последующей фильтрации. Кроме того, это несколько снижает жесткие требования к элементной базе по частотным характеристикам. Но правило остается, чем выше быстродействие аппаратной логики, тем больший диапазон наблюдения можно охватить, не прибегая к различного рода программным ухищрениям и ограничениям при выборе частот преобразования по зонам Найквиста, дискретности, тактовой частоты и даже питания и многое–многое другое.

Промышленные платы Hammerhead от Bittware с успехом справляются с этими условиями. Аппаратную часть виртуального оборудования подключают к промышленному компьютеру, как правило, через шины USB или PCI. Первый вариант не требует вскрытия компьютера, а второй дает обмен на порядок быстрее (спецификация USB3.0 увы пока редко встречается). Задача верхнего уровня сводится к окончательному анализу полученных данных с сочетанием сервисного удобства персонального компьютера (ПК). Кроме того, можно выделить лишь небольшое количество фирм, производящих комплекты для создания виртуального оборудования для работы с DSP в области ВЧ/СВЧ. Наиболее крупные из них – Analog Devices, Bittware, Kontron, National Instruments и Texas Instruments.

Аппаратная часть. Краткое описание объекта

Аппаратной основой (см. рис.2-3) виртуального прибора служит 8–ми слотовое шасси VD3U от компании Kontron [2] с промышленным контроллером CP–306, с установленными периферийными cPCI (Compact PCI) платами Hammerhead фирмы Bittware формата 3U* с четырьмя сигнальными процессорами SHARC ADSP–21160 от Analog Devices. Платы обеспечивают прием и обработку в реальном времени (единицы микросекунд) данных от 4–х каналов АЦП с частотой выборки 32 МГц каждый. Питание обеспечивает встроенный в шасси источник +3.3V/+5V/+12V/–12V. Обмен данными осуществляется по PCI шине. Программная оболочка виртуального прибора в комплексе с аппаратной служит для оценки работоспособности, уровня шумов, джиттера, динамического диапазона модулей АЦП.

* 1U – принятая высота корпуса 44 мм

Предпосылки реализации ПО

Как же получить данные с DSP платы? Все достаточно просто. Для доступа к периферии производитель предоставляет набор драйверов и библиотеку <Hil.dll> или <Hil32v60.dll> в зависимости от версии, подробное описание API которой вы можете узнать в документации [2, 3], а полный список экспортируемых функций вы можете просмотреть в модуле <dsp.pas> (см. ресурсы к статье). Нам-же понадобятся следующие:

• dsp21k_open_by_id()         – массив указателей на область памяти каждого процессора
• dsp21k_reset_bd()              – реинициализация программы в процессоре
• dsp21k_load_exe()             – загрузка прошивки в процессор
• dsp21k_start()                    – старт программы в процессоре
• dsp21k_dl_int()                   – запись значения в область памяти процессора
• dsp21k_ul_int()                   – считывание значения с области памяти процессора
• dsp21k_ul_32s()                 – считывание значения с области памяти процессора (32 бит)

В основу работы программы положен все тот-же алгоритм быстрого преобразования Фурье (БПФ), применяемый к полученным отсчетам сигнала на нижнем уровне и переданными для обработки наверх.

Таким образом, уже можем определить основные требования к нашему виртуальному прибору:

• возможность загрузки прошивки в сигнальные процессоры платы Hammerhead
• визуализация первичных отсчетов, полученных промышленной платой с модуля АЦП
• построение спектра (БПФ/ДПФ) в реальном времени и с возможностью сохранения дампов
• возможность распечатки отображаемых (сохраненных) данных
• управление основными режимами с пульта ДУ

Практика. Разработка ПО и средства отладки

Итак, приступим к основной задаче. Для работы нам понадобится следующее:

• IDE среда разработки Borland Delphi 5–7 (использована для разработки ПО верхнего уровня)
• промышленная плата Hammerhead и модулем АЦП
• генератор сигналов типа Г4–301 или любой другой до 100 МГц
• USB приемник дистанционного управления из материала [4] (см. рис.4)

• любой ИК пульт дистанционного управления для тестирования

Работа с DSP платой

Внешний вид программы реализован через скиновую систему, более подробно (см. исходники). Рассмотрим основные ключевые моменты по доступу к данным… Прежде всего, необходимо подключить и проинициализировать библиотеку для работы с драйвером (см. листинг 1):

ЛИСТИНГ–1

получаем доступ к драйверу

…

function LoadMyDll: bool;

var HLib: THandle;

NewLib: Boolean;

begin

result:= false;

Hlib     := LoadLibrary(’Hil.dll’); // осуществляем динамическое подключение–

NewLib:= (HLib<>0);

if not NewLib then HLib:= LoadLibrary(’Hil32v60.dll’);

if Hlib = 0 then exit;

dsp21k_open_by_id   := GetProcAddress(HLib, ‘dsp21k_open_by_id’);

dsp21k_close        := GetProcAddress(HLib, ‘dsp21k_close’);

if NewLib then begin

dsp21k_load_exe    := GetProcAddress(HLib, ‘dsp21k_load_exe’);

dsp21k_dl_exe      := dsp21k_load_exe;

end else begin

dsp21k_load_exe    := GetProcAddress(HLib, ‘dsp21k_dl_exe’);

dsp21k_dl_exe      := dsp21k_load_exe;

end;

dsp21k_start       := GetProcAddress(HLib, ‘dsp21k_start’);

dsp21k_reset_bd     := GetProcAddress(HLib, ‘dsp21k_reset_bd’);

dsp21k_ul_32s       := GetProcAddress(HLib, ‘dsp21k_ul_32s’);

dsp21k_dl_int       := GetProcAddress(HLib, ‘dsp21k_dl_int’);

dsp21k_ul_int       := GetProcAddress(HLib, ‘dsp21k_ul_int’);

dsp21k_proc_running := GetProcAddress(HLib, ‘dsp21k_proc_running’);

if @dsp21k_open_by_id = nil then exit; // обрабатываем исключения–

if @dsp21k_close = nil then exit;

if @dsp21k_dl_exe = nil then exit;

if @dsp21k_start = nil then exit;

if @dsp21k_reset_bd = nil then exit;

if @dsp21k_dl_32s = nil then exit;

if @dsp21k_dl_int = nil then exit;

if @dsp21k_loaded_file = nil then exit;

if @dsp21k_ul_int = nil then exit;

result:= true

end;

После чего, заведем массивы для хранения отсчетов комплексных и квадратурных составляющих сигнала dim[], qcos[], qsin[] и расчетных значений спектра bpf[]. И инициализируем процессорную плату (см. листинг 2):

ЛИСТИНГ–2

инициализируем процессорную плату

…

Dim      : array [0..10000] of real;

qcos     : array [0..10000] of extended;

qsin     : array [0..10000] of extended;

bpf      : array [0..20000] of extended;

procedure TForm1.init_adc;

var i,j: integer;

p : pointer;

begin

adr_en_data   := $50000;  // адрес флага разрешения передачи для нижнего уровня

adr_val_data  := $52713;  // адрес области памяти с данными

cnt_data      := 1000;   // кол–во отсчетов

id_mod        := 0;      // 0 – слот для платы сбора данных

id_proc       := 1;      // 1 – процессор на плате

try

SetLength(ArrData, cnt_data); // задаем массив для отсчетов

// ===============================================================================

// так как мы заранее не знаем, в каком слоте и в каком процессоре загружена

// прошивка, то “пробегаемся” по всем слотам и процессорам и считываем

// возвращаемый указатель (хэндл)

// в двумерный массив–матрицу indata[][], после чего осуществляем сброс, загрузку

// прошивки и старт процессора на плате Hammerhead

// ===============================================================================

for i:=0  to 63 do //

for j:=1  to 4 do begin

p:= dsp21k_open_by_id(i, j);

if p<>nil then indata.PHandle:= p // просто двумерный

end;

dsp21k_reset_bd(indata[id_mod, id_proc].PHandle);           // сброс прошивки

dsp21k_load_exe(indata[id_mod, id_proc].PHandle, ‘data\21160.dxe’);// загрузка

// ELF файла

dsp21k_start(indata[id_mod, id_proc].PHandle);              // старт прошивки

sleep(10) // делаем на всяк пожарный задержку–

except show_tn(2, ‘Ошибка доступа. Нет модуля или связи…’,’SPEKTRA’); gl_adc:= false end

end;

Получение квадратур сигнала позволяет судить о фазовых неравномерностях в каналах и позволяет производить более тонкую подстройку модулей АЦП. Теория их формирования выходит за рамки данного материала и вряд–ли заинтересует читателя, поэтому приведем лишь сам код (см. листинг 3):

ЛИСТИНГ–3

вычисление квадратур сигнала

…

procedure quadr(auto:boolean;m:integer;var qcos,qsin:array of extended;var dim:array of real;var quad:integer);

var k1, k2, k3, k4, k5, k6, k7, k8, kcos, ksin,

a, aa, b, bb, c, cc, d, dd, f, vcos, vsin: real;

i: integer;

begin

k1:= 0.1169777; // коэффициенты окна–

k2:= 0.4131758;

k3:= 0.7499999;

k4:= 0.9698463;

k5:= 0.9698463;

k6:= 0.7499999;

k7:= 0.4131758;

k8:= 0.1169777;

kcos:= 4*(k1 – k3 + k5 – k7 – k2 + k4 – k6 + k8);

ksin := 4*(k1 – k3 + k5 – k7 + k2 – k4 + k6 – k8);

randomize;

if auto then // формирование псевдопоследовательности

for i:=0 to m do dim:= trunc((512 * sin((sink * i * pi/2) + pi/4) + random(4)) / 4);

m:= ceil(m/32);                // отсчет = 32 бита

quad:= m;                      // кол–во отсчетов

for i:=1 to m do begin

a:= (dim  + dim + dim + dim)*k1;

b:= (dim + dim + dim + dim)*k2;

c:= (dim + dim + dim + dim)*k3;

d:= (dim + dim + dim + dim)*k4;

aa:= (dim + dim + dim + dim)*k1;

bb:= (dim + dim + dim + dim)*k5;

cc:= (dim + dim + dim + dim)*k6;

dd:= (dim + dim + dim + dim)*k7;

f         := a – c + aa – cc;

vcos   := f – b + d – bb + dd;

vsin    := f + b – d + bb – dd;

qcos:= vcos + kcos * 10; // Re – действительная часть

qsin := vsin + ksin * 10; // Im – мнимая часть

end

end;

И собственно то, ради чего все задумывалось, выборка данных с нижнего уровня, передача отсчетов в уже знакомую нам процедуру БПФ, поиск максимума и определение среднего уровня шума в полученном спектре (см. листинг 4):

ЛИСТИНГ–4

передача флага разрешения на нижний уровень и синхронное чтение блока данных

…

procedure TForm1.load_adc;

var a : array [1..20000] of double;

b : array [1..20000] of double;

st: array [0..10000] of real;

exp, mant, lab: string;

i, nf, nn, ns, n: integer;

signoise, re, im, l, m, druc, dmax: extended;

begin

series2.Clear; series3.Clear; inwav1.Clear; // очищаем серии при каждой выборке

dsp21k_dl_int(indata[id_mod, id_proc].PHandle, adr_en_data, 1); //set bit enable

// read data

while dsp21k_ul_int(indata[id_mod, id_proc].PHandle, adr_en_data)<>0  do ;

// ожидаем готовности платы

// передаем на нижний уровень указатель на массив ArrData для заполнения–

dsp21k_ul_32s(indata[id_mod, id_proc].PHandle, adr_val_data, cnt_data,

@ArrData[0]);

for i:= 0 to length(ArrData) – 1 do begin

inwav1.add(ArrData);

dim:= ArrData;    // заносим значения отсчетов в массив dim[] для

// расчета квадратур

end;

if not gl_viborka then begin // выбор режима визуализации – отсчеты / спектр

series2.Assign(inwav1);

series2.SeriesColor:= clred // fix – цвет серии пропадает

end else begin

// ––––––––––––––––––––––––––––

quadr(false, 999, qcos, qsin, dim, quad);

signoise:= –1000; // fix – задаем минимальный уровень шума

nf:= 999;

ch.Title.Text.Text:= ‘Спектр сигнала на выходе приемного модуля при ‘ +

inttostr(nf–1) + ‘–точечном БПФ’;

nn:= ceil(20000000/nf);

for i:=0 to nf do begin

a:= dim;

b:= 0

end;

re:=0; im:=0;

// ==============================================================================

// получение спектра

// ==============================================================================

fft(a, b, nf, 4, 1);

for i:=1  to nf–1 do begin

st:=sqrt(a*a+b*b);  // получаем модуль

if st=0 then st:= 1e–100;

bpf:= 20 * log10(st / nf)    // переводим в дБ

end;

m:=0;

for i:=0 to ceil(nf/2–1) do begin

mant:= inttostr(trunc(i*nn/1000000));

ns    := trunc((frac(i*nn/1000000))*10000);

exp  := inttostr(ns);

if ns<10 then               lab:= mant +’ ,’ + ‘000′ + exp;

if (ns<100)and(ns>9) then   lab:=mant + ‘,’ + ‘00′ + exp;

if (ns<1000)and(ns>99) then lab:=mant + ‘,’ + ‘0′ + exp;

if ns>1000 then             lab:=mant + ‘,’ + exp;

series2.addy(bpf, lab, clblue);

if (series2.YValue > signoise) and (series2.YValue <= 0) then begin

signoise:= series2.YValue;

n:= i

end

end;

// ==============================================================================

// для нормирования по уровню необходимо определить максимум, для этого сканируем

// отсчеты и определяем уровень больший заданного signoise = –1000

// ==============================================================================

signoise:=–1000;

for i:= 10 to nf div 2–5 do begin

if bpf > signoise then begin

signoise:= bpf;

k:= i

end

end;

// ==============================================================================

// поиск среднего шума

// ==============================================================================

nn:= 0;

re := 0;

druc:= –1000;

for i:= 10  to nf div 2–5 do begin

if i < k–10 then begin

re:= re + bpf * bpf;

inc(nn);

if druc < bpf then druc:= bpf

end;

if i > k+10 then begin

re:= re + bpf * bpf;

inc(nn);

if druc < bpf then druc:= bpf

end

end;

re:= –sqrt(re/nn);

// ==============================================================================

// отрисовываем

// ==============================================================================

for i:=0 to ceil(nf/2–1) do series3.Addy(re, ”, clred);

end;

Скомпилировав проект, запустим его на выполнение. Подключив генератор сигналов на вход приемного модуля, подадим тестовый синусоидальный сигнал. В результате уже можем наблюдать отсчеты и при необходимости сами квадратуры (чередующиеся мнимые и действительные составляющие) сигнала (см. рис.5 и 6):

Работа с пультом ДУ

Несмотря на то, что данная «фича» не является основной функцией в такого рода программах, а в некоторых случаях и вредной , но обойти ее стороной никак не могу. Что для этого нужно? Да всего ничего, либо собрать приемник на COM порт и управлять через WinLirc, либо использовать USB.IR приемник, что даст гораздо более стабильные результаты. Было реализовано оба варианта. Чтобы не увеличивать код, сигнатуры нескольких кнопок с ИК пульта были сняты заранее и введены в виде констант. Сама конструкция USB.IR приемника, программа и алгоритм декодирования пакетов были подробно рассмотрены в статье [3], поэтому тут заострять внимание на них не будем и перейдем сразу к коду (см. листинг 5):

ЛИСТИНГ–5

управление режимами с пульта через приемник USB

…

procedure tform1.ic(i: integer);

begin

case i of // нажатия на кнопки

1: image1.OnClick(nil); // чтение с текущего устройства–

2: image2.OnClick(nil); // загрузить файл WAV/MP3 или потоковые данные (дамп)–

3: image3.OnClick(nil); // переключение режимов ”спектр / отсчеты”–

4: image4.OnClick(nil); // вывод на печать–

5: image5.OnClick(nil); // информация о программе–

6: image6.OnClick(nil)  // закрыть программу–

end

end;

// ===============================================================================

// сигнатуры кнопок 1, 2, 3, 4, 5, 6 пульта SR–003

// ===============================================================================

var remote: array[0..5] of string = (

‘12–08–13–07–07–07–13–08–12–08–06–07–06–07–06–08–05–08–06–07–14–07–05–09–05–07–14–06–14–07–12′,

‘13–07–14–07–13–08–05–08–13–07–06–08–13–07–06–07–13–08–06–07–06–07–14–07–06–07–13–08–06–07–12′,

‘13–08–13–07–13–08–13–07–13–07–07–07–06–07–06–08–05–08–06–07–06–07–06–08–06–07–13–08–12–08–11′,

‘13–08–12–08–13–07–13–08–13–08–06–07–13–08–05–08–06–07–06–07–06–07–06–08–05–08–13–07–06–08–14′,

‘13–08–12–07–06–07–06–07–14–08–05–08–05–08–05–08–13–07–06–08–13–07–13–08–06–08–12–07–14–07–12′,

‘14–07–13–07–06–07–06–08–06–07–13–08–13–07–06–08–05–08–06–07–14–07–13–07–13–08–06–07–06–07–12′);

procedure TForm1.tvd3uTimer(Sender: TObject);

var i: integer;

s: string;

j,k,z: smallint;

p: boolean;

begin

if (n5.Checked)and(gl_adc) then load_adc;  // визуализация спектра в динамике–

// ===============================================================================

// получение кода кнопки с USB приемника (сигнатуры)

// ===============================================================================

if (GetInfraCode(’ra_usb’) = 1) then Exit;

if DataLength = 0 then Exit;

for i:= 0 to DataLength–1 do begin

s:= s + inttohex(InputInfraData, 2);

if i<> DataLength–1 then s:= s + ‘–’

end;

// ===============================================================================

// проверка сигнатуры с учетом интервала доверия и передача на

// интерпретатор команд

// ===============================================================================

for i:=1 to 6 do

for z:= 0 to DataLength–1 do

if (strtoint(’$’+copy(s,(z*3)+1,2))–dover >=

strtoint(’$’+copy(remote,(z*3)+1,2)))or

(strtoint(’$’+copy(remote,(z*3)+1,2)) <=

strtoint(’$’+copy(s,(z*3)+1,2))+dover) then begin

p:= true;

k:= z

end else begin p:= false; break end;

if p then ic(k)

end;

Осталось проверить работоспособность управления по ИК. Для этого, подключив USB.IR приемник, понажимаем кнопку «2» на пульте, отвечающую за смену режима отображения «спектр / отсчеты» (см. рис.7):

Заключение

Рассмотренный виртуальный прибор позволяет сэкономить время и снизить стоимость любой системы сбора и анализа данных с модулей АЦП без привлечения дорогой специализированной аппаратуры, за счет применения гибкого ПО в сочетании с производительностью DSP и простоты «Plug&Play» подключения (cPCI) промышленных плат фирмы Bittware.

Полные исходные тексты и компиляцию виртуального спектроанализатора SPEKTRA (файл fft2.zip) вы можете загрузить на форуме клуба программистов (раздел «Журнал клуба программистов.  Первый выпуск») или с сайта автора [5]. Если тема представляет для вас интерес – пишите, задавайте вопросы на форуме  http://www.programmersforum.ru

Ресурсы

• С.Бадло. Быстрое преобразование Фурье. Практика использования. – Блог клуба программистов,

05.02.2010 http://pblog.ru/?p=658 и http://www.programmersforum.ru/showthread.php?t=83467

• Data Sheet Kontron Modular Computers GmbH, 2003, ID26799, rev.01
• Practical Design Techniques for Sensor Signal Conditioning, Analog Devices, 1998
• E.Бадло, С.Бадло. USB термометр и дистанционка в одном флаконе. Часть 2. – Радиолюбитель, 2010,

№1, с.48 http://raxp.radioliga.com/cnt/s.php?p=us2.pdf или с форума клуба программистов

http://www.programmersforum.ru/attachment.php?attachmentid=17684&d=1258320468

• Ресурсы и компиляция проекта SPEKTRA http://raxp.radioliga.com/cnt/s.php?p=fft2.zip

Статья из первого выпуска “журнала ПРОграммистов”.
Скачать этот номер можно по ссылке.
Ознакомиться со всеми номерами журнала.

Обсудить на форуме — Быстрое преобразование Фурье. Практика использования. Часть 2

Похожие статьи

• Компилятор домашнего приготовления. Часть 4
• Общение между запущенными копиями своих программ
• Близкие контакты третьего вида с Visual Foxpro
• Запись дисков в Delphi средствами IMAPI
• WMI. Wладение Mагической Iнформацией. Часть 2
• Быстрое написание программ на WinAPI
• БПФ. Практика использования
• WMI. Wладение Mагической Iнформацией. часть 1
• Программы для самых маленьких
• Делаем динамические тени на OPENGL. Часть 1