Последние записи

Интенсив по Python: Работа с API и фреймворками 24-26 ИЮНЯ 2022. Знаете Python, но хотите расширить свои навыки?
Slurm подготовили для вас особенный продукт! Оставить заявку по ссылке - https://slurm.club/3MeqNEk

Online-курс Java с оплатой после трудоустройства. Каждый выпускник получает предложение о работе
И зарплату на 30% выше ожидаемой, подробнее на сайте академии, ссылка - ttps://clck.ru/fCrQw

23rd
Фев

БПФ. Практика использования

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

Получение спектра в радиотехнике уже стало обыденным явлением. Появились как аппаратные высокоскоростные реализации, например от таких брендов как Tektronix, так и совмещенные варианты анализаторов на основе DSP процессоров или ПЛИС в промышленных или офисных компьютерах. Данным материалом мы начинаем цикл статей посвященных теме анализа спектра сигналов и их визуализации, для чего сегодня разработаем компонент, работающий с цифровым аудиопотоком, и освоим методику Фурье-анализа применительно к распознаванию DTMF.

Сергей Бадло
by raxp http://raxp.radioliga.com

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

Спектроанализатор* — это прибор для наблюдения и измерения относительного распределения энергии электромагнитных колебаний в заданной полосе частот и бывает как параллельного или последовательного типа, так и совмещенным. По способу обработки — различают аналоговые и цифровые, а по характеру анализа — скалярные (получение частотно-амплитудных спектров) и векторные (фазо-частотных спектров).

Анализатор спектра позволяет определить амплитуду и частоту спектральных составляющих, входящих в состав анализируемого сигнала. Важнейшим его параметром – является разрешающая способность, т.е. наименьший интервал по частоте между двумя гармониками, которые еще можно измерить.

Рис. 1. «Преимущества софтовых вариантов очевидны лишь на малых частотах, либо при использовании аппаратно-программных реализаций»

Физический смысл или… для чего мы учим математику

Вспомним курс математики [1-6]. Как вы знаете, периодическим сигналом называют такой вид воздействия, когда форма сигнала повторяется через некоторый интервал времени T, который называется периодом. Простейшей формой периодического сигнала является гармонический сигнал или синусоида, которая характеризуется амплитудой, фазой и периодом. Все остальные сигналы будут негармоническими.

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

U(t) = Uo + SUM ( Um * sin ( k * ? * t + ? ) );

где k — номер гармоники;
k? — угловая частота k- ой гармоники;
? = 2*pi/T — угловая частота первой гармоники;
? — начальная фаза сигнала;
Uo — нулевая гармоника.

Для выделения спектра в радиотехнике, как правило, используется быстрое преобразование Фурье (БПФ). БПФ — это быстрый алгоритм вычисления дискретного преобразования Фурье. То есть алгоритм вычисления за количество действий, меньшее чем O(N2), требуемых для прямого вычисления ДПФ.

Для чего нужно быстрое преобразование Фурье? Допустим у нас есть периодическая функция изменяющаяся по закону синуса x = sin(t) (см. рисунок 2). Максимальная амплитуда этого колебания равна 1. Если умножить его на некоторый коэффициент A, то получим тот же график (см. риснок 3), растянутый по вертикали в A раз: x = A*sin(t)

Рис. 2. Периодическая функция

Рис. 3. Увеличение амплитуды

Период колебания равен 2pi. Если мы хотим увеличить период до T, то надо умножить переменную t на коэффициент [0; 1]. Это вызовет растяжение графика по горизонтали: x = A*sin(2pi/T). Как вы знаете, частота колебания обратна периоду: f = 1/T. Также говорят о круговой частоте, которая вычисляется по формуле: ? = 2pi/T, где x = A*sin(?t).

И, наконец, есть фаза, обозначаемая как ?. Она определяет сдвиг графика колебания влево. В результате сочетания всех этих параметров получается гармоническое колебание (гармоника) или спектральная составляющая.

Если изменить фазу на 90 градусов, то можно перейти от синуса к косинусу. Для удобства, далее будем работать с функцией косинуса:

x = A * cos ( 2pi/T + ? ) = A * cos ( ?t + ? );

Преобразуем по формуле косинуса суммы:

x = A * cos ? * cos ( 2pi/T ) — A * sin ? * sin ( 2pi*t/T );

Выделим элементы, независимые от времени t, и обозначим их как Re и Im (действительная и мнимая части):

x = Re * cos ( 2pi*t/T ) – Im * sin ( 2pi*t / T );
Re = A * cos ?, Im = A * sin ?;

По величинам Re и Im можно однозначно восстановить амплитуду и фазу исходной гармоники:

? = arctg ( Im/Re ), A = sqrt ( Re^2 + Im^2 );

Теперь возьмем обратное преобразование Фурье:

Xn = ( 1/N ) * SUM ( Xk* e ^ ( j*2*pi*k*n / N ) );

И выполним над этой формулой следующие действия:

разложим каждое комплексное Xn на мнимую и действительную составляющие Xn = Re + j*Im разложим экспоненту по формуле Эйлера на синус и косинус действительного аргумента перемножим внесем множитель 1/N под знак суммы и перегруппируем элементы в две суммы:

Xn = ( 1/N ) * SUM ( Xk * [ cos ( 2*pi*k*n/N ) + j * sin ( 2*pi*k*n/N ) ] ) =>
=> ( 1/N ) * SUM ( ( Rek + j*Imk ) * [ cos ( 2*pi*k*n/N ) + j*sin ( 2*pi*k*n/N ) ] ) =>
=> SUM ( ( Rek/N )* cos ( 2*pi*k*n/N ) – ( Imk/N ) * sin ( 2*pi*k*n/N ) ) +
+ j*SUM ( ( Rek/N )* sin ( 2*pi*k*n/N ) + ( Imk/N ) * cos ( 2*pi*k*n/N ) );

Как видите, слева стоит действительное число Xn, а справа две суммы, одна из которых помножена на мнимую единицу j. Сами же суммы состоят из действительных слагаемых. Отсюда следует, что вторая сумма равна нулю, если исходная последовательность была действительной. Отбросим ее и получим:

Xn = SUM ( ( Rek/N ) * cos ( 2*pi*k*n/N ) – ( Imk/N ) * sin ( 2*pi*k*n/N ) );

Поскольку при дискретизации мы брали tn = nT/N и Xn = F(tn), то можем выполнить замену: n = tn*N/T. Следовательно, в синусе и косинусе вместо 2pi*k*n/N можно написать 2pi*k*tn/T. В результате получим:

Xn = F(tn) = SUM ( ( Rek/N ) * cos ( 2*pi*k*tn/N ) – ( Imk/N ) * cos ( 2*pi*k*tn/N ) );

Сопоставим эту формулу с формулами для гармоники:

x = A * cos ( 2*pi*t/T + ? ) = A * cos ( ?t + ? );
x = Re * cos ( 2*pi*t/T ) – Im * sin ( 2pi*t / T );

Следовательно, сумма представляет собой сумму из N гармонических колебаний разной частоты, фазы и амплитуды:

F(tn) = SUM ( Ak * cos ( 2*pi*tn/Tk + ?k ) ) = SUM ( Gk(tn) );

Далее будем функцию Gk(t) = Ak*cos(2pi*tk/T + ?k) называть k-й гармоникой. Амплитуда, фаза, частота и период каждой из гармоник связаны с коэффициентами Xk формулами:

Xk = Rek + j * Imk;
Xk = N * Ak * e ^ ( j * ?k );
Ak = ( 1/N ) * sqrt ( Rek^2 + Imk^2 );
?k = arctg ( Imk / Rek );
Tk = T/k;

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

Описание компонента спектроанализатора

Разрабатываемый компонент предназначен для построения спектра аудио-сигнала, кодирования и декодирования двух-тоновых посылок DTMF (Dual Tone Multi Frequency) и получения «сырых» отсчетов в реальном времени. Его можно использовать в системах сигнализации, различных плеерах аудио-видео файлов и учебных программах работы со звуком. В основу работы компонента положено использование алгоритма быстрого преобразования Фурье (БПФ).

Входным является внутренний буфер с аудиоданными, частотой дискретизации 44100 герц и форматом 16 бит/семпл. Длина буфера фиксирована, в данной версии компонента выбор не реализован и ограничен величиной в 3000 отсчетов. Сам компонент невизуальный.

Внешние свойства и события компонента:

property About — Copyright 🙂
property DTMF_keys — строка для генерации DTMF
property DTMF_volume — установка амплитуды генерации
property DTMF_duration_ms — установка длительности генерации
property FFT_point — выбор количества точек преобразования БПФ
property FFT_window — выбора типа сглаживающих окон
property Key — событие декодированных команд DTMF
property Spektra — // — выдача спектра после БПФ
property DataOsc — // — выдача «сырых» отсчетов с аудио-буфера

Результат работы компонента и типичный спектр сигнала DTMF с его распознаванием представлен на (см. рисунок 4):

Рис. 4. Визуализация сигнала. Типичный спектр сигнала DTMF

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

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

среда разработки TurboDelphi-Lite portable
аудиокарта

Вкратце, процедура (прямого) БПФ в компоненте будет включать в себя следующие шаги:

берем из сигнала N выборок кратным степени 2, т.е. 2^k
рассчитываем комплексное БПФ, мнимые части заполняем нулями, получаем 2N значений
амплитуду сигнала для каждой гармоники получаем складывая квадраты действительной и мнимой части и извлекая из суммы корень квадратный
получаем N значений, из которых значения от 0 до (N/2-1) представляют наш спектр в области от 0 до половины частоты дискретизации, вторую половину (зеркалку) отбрасываем
для адекватного представления пересчитываем в дБ, с учетом максимальной величины в выборке по формуле 20lg(Ai/Amax), для напряжений
при необходимости используем различные сглаживающие окна для взвешивания входного сигнала во временной области, например Блэкмана-Харриса
добавляем порог чувствительности (подставку)
результаты выводим в качестве события компонента, например, используя series для подключения к TChart-у

Разбор принципов генерации и декодирования DTMF сигналов проведен в статье [6] и в данном материале рассматриваться не будет. В листинге 1 приведен полный код компонента с подробными комментариями…

…
unit DTMF;
interface
uses MMSystem, Windows, SysUtils, Messages, Classes, controls, extctrls, series, TeEngine, math;
type// тип данных wave- ind
  TData16 = array[0..127]ofsmallint;
  PData16 = ^TData16;
Type// установки для waveform
 SINEWAVE = packedrecord
  dblFrequency : Double;
  dblDataSlice : Double;
  dblAmplitudeL : Double;
  dblVolumeF : Double;
end;
type
  Twindow = (dB_0,
                    dB_54,
                    dB_67,
                    dB_72,
                    dB_92); // функции окна-
type
  Tkeys = procedure(Sender:TObject; key:string; a1,a2,f1,f2: double)ofobject; // выдача DTMF
  Tspektr = procedure(Sender:TObject; series: TbarSeries)ofobject; // выдача спектра
  TdataOsc = procedure(Sender:TObject; series: TfastlineSeries)ofobject; // сырой набор данных
  TDTMF=class(TComponent)
private
    fabout : string;
    fkey : string;
    fvol : integer;
    flen : integer;
    fcntp : integer;
    FOnKeys : TKeys;
    FOnSpektr : TSpektr;
    FOnDataOsc : TDataOsc;
    ftimer : Ttimer;
    ftimer2 : Ttimer;
    tmr_en : boolean;
    fwindow : twindow;
protected
procedure gen_dtmf(const Value: string); // передача строки DTMF для генерации
procedure setabout(const Value: string); // мой Copyright 🙂
procedure set_window(const Value: twindow); // выбор окна сглаживания
procedure f_cntp(const Value: integer); // установка к-ва точек БПФ
procedure wcard; // инит-деинит работы с аудио
procedure ind(Sender: TObject); // события компонента
procedure ind2(Sender: TObject); // генерация одиночного DTMF с перебором
public
constructor Create(AOwner: TComponent); override;
destructor Destroy; override;
published// внешние свойства компонента
property About : stringread Fabout write setabout;
property DTMF_keys : stringread Fkey write gen_dtmf;
property DTMF_volume : integerread Fvol write fvol default 100;
property DTMF_duration_ms : integerread Flen write flen default 250;
property FFT_point : integerread Fcntp write f_cntp default 2048;
property FFT_window : twindow read fwindow write set_window;
property Key : TKeys read FOnKeys write FOnKeys;
property Spektra : TSpektr read FOnSpektr write FOnSpektr;
property DataOsc : TDataOsc read FOnDataOsc write FOnDataOsc;
end;
procedureRegister;
const// таблица соответствия частот DTMF
    keys = ‘1234567890*#abcd’;
    DTMF1: array[1..16]ofinteger
=(697,697,697,770,770,770,852,852,852,941,941,941,697,770,852,941);
    DTMF2: array[1..16]ofinteger
=(1209,1336,1477,1209,1336,1477,1209,1336,1477,1336,1209,1477,1633,1633,1633,1633);
var stp: boolean = FALSE;
    inwav, outwav : TfastLineSeries;
    spektr : TbarSeries;
// декодер DTMF
    adr2 : pWaveHdr;
    BufHead1,BufHead2 : TWaveHdr;
    bufsize : integer;
    header : TWaveFormatEx;
    hwi2 : HWAVEIN;
    hBuf : THandle;
    pnt : PPointArr;
    f_window : smallint = 5; // тип окна – без сглаживания
    fcntpp : integer = 2048; // к-во точек FFT
    signal : string; // декодированный DTMF (временная переменная)
    a1,a2, // амплитуды-
    f1,f2 : double; // частоты-
    // кодер DTMF
    waveOut : hWaveOut;
    outHdr : TWaveHdr;
    header2 : TWaveFormatEx;
    pBuf : tHandle;
    pBuffer : pointer;
    Opened, lock : boolean;
    gl_key : integer;
implementationSyhi-подсветка кода

Обычно, генерацию звука в памяти и воспроизведение в среде Windows осуществляют через Waveform Audio Win32 API. Нам понадобятся следующие функции:

waveOutOpen — открывает имеющееся устройство вывода Waveform audio для сигнала
waveOutPrepareHeader — выполняет подготовку буфера для операции вывода данных

Далее зададимся законом модуляции, форматом вывода PCM, частотой дискретизации, количеством каналов и длительностью генерации…

//——————————————————————————
// КОДЕР DTMF – генерация 2-х тонального сигнала
//——————————————————————————
procedure stopPlay;
begin
 WaveOutReset(WaveOut);
 WaveOutClose(WaveOut);
 GlobalUnlock(pBuf);
 GlobalFree(pBuf)
end;
procedure PlayBuffer(h: hwnd; SampleRate: integer; Bits: Byte; Buffer: array of byte);
var Err: integer;
begin
 with header2 do begin
  wFormatTag := WAVE_FORMAT_PCM;
  nChannels := 1;
  nSamplesPerSec := SampleRate;
  wBitsPerSample := Bits;
  nBlockAlign := nChannels * (wBitsPerSample div 8);
  nAvgBytesPerSec := nSamplesPerSec * nBlockAlign;
  cbSize := 0;
 end;
 if Opened = true then stopPlay;
 err:= WaveOutOpen(addr(waveOut), 0, @header2,h, 0, CALLBACK_WINDOW);
 if Err <> 0 then Exit;
 pBuf := GlobalAlloc(GMEM_MOVEABLE and GMEM_SHARE, length(Buffer));
 pBuffer := GlobalLock(pBuf);
 with outHdr do begin
  lpData := PBuffer;
  dwBufferLength := length(Buffer);
  dwUser := 0;
  dwFlags := 0;
  dwLoops := 0
 end;
 err:= WaveOutPrepareHeader(waveOut, @outHdr, sizeof(outHdr));
 if Err <> 0 then Exit;
 copyMemory(pBuffer, @Buffer, length(Buffer));
 err:= WaveOutWrite(waveOut, @outHdr, sizeof(outHdr));
 if Err <> 0 then Exit
end;
function sel_byte(lngWord: Longint; intPosition: Smallint): byte;
var lngTemp: Longint;
    intByte: byte;
begin
 if intPosition=3 then begin
 // Byte 2
 lngTemp := lngWord;
  // Mask off byte and shift right 24 bits.
  // Mask -> 2130706432 = &H7F000000
  // Shift -> Divide by 16777216
  lngTemp := Round((lngTemp and 2130706432)/16777216);
  intByte := lngTemp;
 end else if intPosition=2 then begin
  // Byte 2
  lngTemp := lngWord;
  lngTemp := Round((lngTemp and 16711680)/65536);
  intByte := lngTemp;
  end else if intPosition=1 then begin
   // Byte 1
   lngTemp := lngWord;
  // Mask off high byte and shift right 8 bits.
  // Mask -> 65290 = &HFF00
  // Shift -> Divide by 256
  lngTemp := Round((lngTemp and 65290)/256);
  intByte := lngTemp;
  end else begin
   // Byte 0
   intByte := lngWord and $FF;
  end;
 result:= intByte
end;
procedure toneGenerate(lngSampleRate: integer; intBits: byte; dblVolume: array of double; var Freq:
 array of Smallint; Seconds: Double; var FreqBuffer: variant); // создание WAVEFORM
var i, j : integer;
    lngLimit, lngData : longint;
    lngSamples, lngDataSize : integer;
    dblDataPtL, dblWaveTime,
dblSampleTime, dblFrequency: Double;
tmpBuf : array of byte;
intSineCount : Smallint;
SineWaves : array of SINEWAVE;
begin
setLength(SineWaves, length(freq));
for i:=0 to length(freq) — 1 do begin
with SineWaves[i] do begin
dblAmplitudeL:= 0.25;
dblFrequency := freq[i]; // задаем частоты генерации WAVEFORM
dblVolumeF := dblVolume[i]
end
end;
intSineCount := length(SineWaves)-1;
for i:=0 to intSineCount do begin
dblWaveTime := 1 / SineWaves[i].dblFrequency;
dblSampleTime := 1 / lngSampleRate;
SineWaves[i].dblDataSlice := (2*Pi)/(dblWaveTime/dblSampleTime);
end;
lngSamples := round(Seconds/dblSampleTime);
lngDataSize := Round(lngSamples*(intBits/8));
SetLength(tmpBuf, lngDataSize);
if intBits=8 then lngLimit := 127
else lngLimit := 32767;
for i:=0 to lngSamples-1 do begin
if intBits=8 then begin
// ————————————————————————
// 8 Bit Data
// ————————————————————————
dblDataPtL := 0;
for j:=0 to intSineCount do
dblDataPtL := dblDataPtL +
(sin(i*SineWaves[j].dblDataSlice)*SineWaves[j].dblAmplitudeL*SineWaves[j].dblVolumeF);
lngData := round(dblDataPtL*lngLimit)+lngLimit;
tmpBuf[i] := ExtractByte(lngData, 0);
end else begin
// ————————————————————————
// 16 Bit Data
// ————————————————————————
dblDataPtL := 0;
for j:=0 to intSineCount do
dblDataPtL := dblDataPtL +
(sin(i*SineWaves[j].dblDataSlice)*SineWaves[j].dblAmplitudeL*SineWaves[j].dblVolumeF);
lngData := Round(dblDataPtL*lngLimit);
tmpbuf[2*i] := sel_byte(lngData, 0);
tmpbuf[(2*i)+1] := sel_byte (lngData, 1);
end
end;
FreqBuffer:= tmpBuf
end;
procedure tdtmf.gen_dtmf(const Value: string); // передача строки DTMF для генерации-
begin
fkey:= value;
if fkey<>» then tmr_en:= true // запрещаем генерацию, если пусто
end;
procedure TDTMF.ind2(Sender: TObject); // генерация одиночного DTMF с перебором-
var Freq: array [0..1] of smallint;
Buffer:array of byte;
dblVolume: array [0..1] of double;
SoundBuffer: variant;
i: integer;
begin
if tmr_en then begin
inc(gl_key);
if gl_key > length(fkey) then begin
gl_key:= 0;
tmr_en:= false; // если перебрали все введенные символы – останов генерации
fkey:= »
end else for i:= 1 to length(keys) do
if keys[i]= lowercase(fkey[gl_key]) then begin
Freq[0]:= dtmf1[i]; // задаем частоты-
Freq[1]:= dtmf2[i];
dblVolume[0]:= fvol / 100; // задаем уровень громкости-
dblVolume[1]:= fvol / 100;
toneGenerate(22050, 8, dblVolume, Freq, flen/1000, SoundBuffer); // создание WAVEFORM-
buffer:= SoundBuffer;
PlayBuffer(0,22050, 8, Buffer) // воспроизведение-
end
end
end;
// END DTMF KODER ————————————————————-Syhi-подсветка кода

Как быть с обработкой в реальном времени? Воспользуемся API функцией WaveInOpen, чтобы получить доступ к текущему аудиоустройству. Также заведем два буфера BufHead1 и BufHead2, один для накопления, второй для получения данных. Размер буфера определим в 3000 отсчетов, т.к. нам не требуется высокое разрешение по частоте, допустимую погрешность при определении DTMF будем задавать доверительным интервалом по частоте. Частоту дискретизации зададим типичную (максимальную) для большинства аудиокарт в 44100 Гц, 16 бит на канал. После чего, передадим полученный набор данных в нашу процедуру БПФ и строим спектр как обычно. Причем, заметьте, основное время тратится не на обработку данных и БПФ, а на набивку в series. Поэтому, если вам дорого время и вы хотите максимально увеличить размер буфера для повышения разрешения по частоте, то придется отказаться от удобства использования TChart (именно этим обусловлено использование series)…

//————————————————————————————
// ДЕКОДЕР DTMF + СПЕКТРОАНАЛИЗ
//————————————————————————————
function FFT(var x, y:arrayofdouble;var nn:integer;nf, ii: integer): integer;
var c,s,t1,t2,t3,t4,u1,u2,u3,z,a0,a1,a2,a3,w:double; // функции окна-
    i,j,p,rt,l,ll,m,m1,k:integer;
begin
 rt:= 1;
 nn:= nn div2;
while rt<nn do
  rt:=rt*2;
 nn:= rt;
 z:= 2*pi/nn;
// выбор окна подавления
…
a0:=1; // без изменений
a1:=0;
a2:=0;
a3:=0;
for i:=0to nn-1dobegin
  w:=a0-a1*cos(z*i)+a2*cos(z*2*i)+a3*cos(z*3*i);
  x[i]:=x[i]*w;
  y[i]:=y[i]*w;
end;
//——————————————-
 ll:= nn;
 M := nn div2;
 M1:= Nn-1;
while ll>=2dobegin
  l:=ll div2;
  u1:=1;
  u2:=0;
  t1:=PI/l;
  c:=cos(t1);
  s:=(-1)*ii*sin(t1);
for j:=0to l-1do
begin
    i:=j;
while i<nn do
begin
      p:=i+l;
      t1:=x[i]+x[p];
      t2:=y[i]+y[p];
      t3:=x[i]-x[p];
      t4:=y[i]-y[p];
      x[p]:=t3*u1-t4*u2;
      y[p]:=t4*u1+t3*u2;
      x[i]:=t1;
      y[i]:=t2;
      i:=i+ll;
end;
    u3:=u1*c-u2*s;
    u2:=u2*c+u1*s;
    u1:=u3;
end;
  ll:=ll div2
end;
 j:=0;
for i:=0to m1-1dobegin
if i>j thenbegin
   t1:=x[j];
   t2:=y[j];
   x[j]:=x[i];
   y[j]:=y[i];
   x[i]:=t1;
   y[i]:=t2;
end;
 k:=m;
while j>=k dobegin
  j:=j-k;
  k:=k div2;
end;
 j:=j+k;
end
end;
procedure FFTQuad(seriesin,seriesout: TChartSeries; max:integer); // max- точечное БПФ
var a,b : array of double;
      i,k : integer;
      d : real;
begin
 i:=0;
 if Seriesin.yValues.count = 0 then exit;
  k:= Seriesin.YValues.Count;
 while (k>1) and (power(2, i)<max) do
  begin
   k:=k div 2;
   inc(i)
  end;
 k:= ceil(power(2, i));
 SetLength(a,k); // инициализируем массив Re, Im
 SetLength(b,k);
 for i:=0 to k-1 do
 begin
  a[i]:= Seriesin.YValue[i];
  b[i]:= 0
 end;
 FFT(a, b, k, f_window, 1); // домножаем на выбранное окно, получение спектра
 for i:=0 to k div 2-1 do begin // отсекаем зеркалку
  d:= sqrt(a[i]*a[i] + b[i]*b[i]); // получаем модуль из Re и Im
  d:= 20*log10(d/k + 0.000001) -25; // приведение к дБ и нормирование
// -25дб это подставка, чтоб убрать фоновые шумы
// по спецификации ITU-T для DTMF
// уровень шума (SNR) на уровне 15 дБ
  seriesout.Add(d)
 end
end;
// получение аудиоданных и построение спектра —————————————————-
procedure waveInProc2(hwi: HWAVEIN; uMsg,dwInstance,dwParam1,dwParam2: DWORD);stdcall;
var i : integer;
      data16 : PData16;
      temp : pWaveHdr;
      a,f,cntval : double;
begin
 if (uMsg=WIM_DATA) and (stp) then begin
  temp:= adr2;
  if adr2= @bufhead1 then adr2:= @bufhead2 // получаем указатель на данные с буфера 1/2
   else adr2:= @bufhead1;
  //
  if stp then WaveInAddBuffer(hwi,adr2,SizeOf(TWaveHdr));
  data16:= PData16(temp.lpData); // собственно сами данные
  if (not lock) then try inwav.Clear; outwav.Clear; spektr.Clear; // подчищаем
  for i := 0 to BufSize — 1 do begin // набиваем из аудиобуфера-
   inwav.add(data16^[i])
  end;
  // ПОЛУЧЕНИЕ СПЕКТРА —
  FFTQuad(inwav, outwav, fcntpp);
  // обработка спектра и 2-x проходный поиск —
  a1:= -1000;
  cntval:= header.nSamplesPerSec / outwav.YValues.Count;
  for i:= 0 to (outwav.YValues.Count)-1 do begin
   a:= outwav.YValues[i];
   f:= i * cntval; // получение истинной частоты гармоники
   if a>=0 then spektr.AddXY(f,a)
    else spektr.AddXY(f,0); // отсекаем отрицательные амплитуды
   if a > a1 then begin a1:= a; f1:= f end // частота для max 1- гармоники
  end;
  a2:= -1000;
  for i:= (outwav.YValues.Count)-1 downto 0 do begin
   a:= outwav.YValues[i];
   f:= i * cntval;
   if (a > a2) and (a<>a1) then begin a2:= a; f2:= f end // частота для max 2- гармоники
  end;
  // ———————————————————————————————————
  // ИДЕНТИФИКАЦИЯ DTMF
// ———————————————————————————————————
// по спецификации ITU-T на DTMF доверительный интервал должен быть не более 1.5%,
// но мы зададимся чуть больше, чтобы учесть разброс характеристик и небольшое заданное
// разрешение анализатора спектра по частоте (размер буфера 3000, см. выше)
signal:= »;
for i:= 1 to 16 do begin
if (DTMF2[i]*0.96<f1) and (DTMF2[i]*1.04>f1) and // 1 амплитуда >2
(DTMF1[i]*0.96<f2) and (DTMF1[i]*1.04>f2) then begin
signal:= keys[i];
break
end;
if (DTMF1[i]*0.96<f1) and (DTMF1[i]*1.04>f1) and // 1 амплитуда >2
(DTMF2[i]*0.96<f2) and (DTMF2[i]*1.04>f2) then begin
signal:= keys[i];
break
end;
end;
spektr.Title:= ‘DTMF(‘+ signal +‘): ‘ +
format(‘A1= %.2n’,[a1])+ formatfloat(‘ [0 Hz] ‘, f1) +
format(‘A2= %.2n’,[a2])+ formatfloat(‘ [0 Hz]’, f2);
except end
end else Exit
end;
//————————————————————————————————————
// инициализация-деинициализация получения аудио-данных
procedure TDTMF.wcard;
const rbuf = 6;
var BufLen : word;
buf : pointer;
begin
stp:= not stp;
try
if stp then begin // старт
BufSize:= rbuf *500;
with header do begin
wFormatTag:= WAVE_FORMAT_PCM;
nChannels := 2; // каналов
nSamplesPerSec:= 44100; // дискретизация, Гц
wBitsPerSample:= 16; // бит
nBlockAlign:= nChannels * (wBitsPerSample div 8);
nAvgBytesPerSec:= nSamplesPerSec * nBlockAlign;
cbSize:= 0;
end;
WaveInOpen(Addr(hwi2), WAVE_MAPPER, addr(header),integer(@waveInProc2),
0,CALLBACK_FUNCTION);
BufLen:= header.nBlockAlign * BufSize;
hBuf:= GlobalAlloc(GMEM_MOVEABLE and GMEM_SHARE, BufLen);
Buf:= GlobalLock(hBuf);
with BufHead1 do begin
lpData:= Buf;
dwBufferLength:= BufLen;
dwFlags:= 0;
end;
with BufHead2 do begin
lpData:= Buf;
dwBufferLength:= BufLen;
dwFlags:= 0;
end;
adr2:= @BufHead1;
waveInPrepareHeader(hwi2, Addr(BufHead1), sizeof(BufHead1));
waveInPrepareHeader(hwi2, Addr(BufHead2), sizeof(BufHead2));
WaveInAddBuffer(hwi2, addr(BufHead1), sizeof(BufHead1));
WaveInStart(hwi2) 
end else begin // стоп
WaveInReset(hwi2);
WaveInUnPrepareHeader(hwi2, addr(BufHead1), sizeof(BufHead1));
WaveInClose(hwi2);
GlobalUnlock(hBuf); GlobalFree(hBuf);
end
//
except end
end;
//————————————————————————————————————
// СОБЫТИЯ КОМПОНЕНТА
//————————————————————————————————————
// в принципе содержимое можно перенести в процедуру waveInProc2(), но хотелось гибкости
//————————————————————————————————————
procedure TDTMF.ind(Sender: TObject);
begin
lock:= true; // блокируем очистку, БПФ и декодирование DTMF на время выдачи
inwav.SeriesColor := rgb(0,0,255); // синий цвет серии
spektr.SeriesColor:= rgb(255,0,0); // красный-
if Assigned(FOnSpektr) then FOnSpektr(self, spektr); // спектр
if Assigned(FOnDataOsc) then FOnDataOsc(self,inwav); // сырые данные (осциллограф)
if signal<>» then
if Assigned(FOnKeys) then FOnKeys(self, signal, a1, a2, f1, f2); // декодированный DTMF
lock:= false
end;
//————————————————————————————————————
// СЕРВИС-МОДУЛЬ (СКЕЛЕТ)
//————————————————————————————————————
constructor TDTMF.Create(AOwner: TComponent);
begin
inherited Create(aowner);
fvol := 100;
flen := 250;
Fcntp := fcntpp; // устанавливаем по умолчанию 2048 точек БПФ
fwindow:= dB_0;
Fabout := ‘Badlo Sergey’;
Inwav := tfastlineseries.Create(nil);
Outwav := tfastlineseries.Create(nil);
Spektr := tbarseries.Create(nil); // визуально удобнее-
spektr.Marks.Visible:= false;
ftimer:= ttimer.Create(self);
ftimer.Enabled := false;
ftimer.interval := 200;
ftimer.ontimer := ind;
ftimer.Enabled := true;
//
ftimer2:= ttimer.Create(self);
ftimer2.Enabled := false;
ftimer2.interval := 350;
ftimer2.ontimer := ind2;
ftimer2.Enabled := true;
wcard // инициализация получения аудиоданных
end;
destructor TDTMF.Destroy;
begin
stopplay; // запрещаем генерацию DTMF
wcard; // деинициализация аудио
ftimer.Free; ftimer2.Free;
inwav.Free; outwav.Free; spektr.Free;
inherited
end;
procedure TDTMF.f_cntp(const Value: integer);
begin
fcntp:= value; fcntpp:= value
end;
procedure TDTMF.set_window(const Value: twindow); // выбираем окно сглаживания-
begin
fwindow:= value;
case value of
dB_0 : f_window:= 5; // без сглаживания
…
end
end;
procedure TDTMF.setabout(const Value: string);
begin
fabout:= ‘Badlo Sergey’
end;
//============================================================
procedure Register;
begin
RegisterComponents(‘RAMEDIA’, [TDTMF])
end;
end.Syhi-подсветка кода

Осталось проверить работоспособность нашего модуля. Для этого (см. рисунки 5-7):

Рис. 5. Инсталлируем компонент

Выбрав в меню «Component/Install Component», проинсталлируем компонент DTMF (модуль <dtmf.pas>). Создав новый проект «File/New Application», перенесем наш компонент на форму. После чего, станут доступными его свойства и методы (см. рис.6). Также нам понадобится TChart (визуализация), TMemo (декодированные посылки) и TEdit (строка ввода команд генерации) (см. рисунок 7). Теперь напишем следующий код в событиях компонента (см. листинг 2):

Рис. 6. Задаем параметры

Рис. 7. События компонента

тестовый модуль
…
procedure TForm1.DTMF1Key(Sender: TObject; key: String; a1, a2, f1,
  f2: Double);
var s: string;
begin
 CAPTION:= ‘ Тестовый спектроанализатор. Декодер-кодер DTMF (‘ + key + ‘)’;
 s:= key +
     format(‘ — A1= %.2n’,[a1])+ formatfloat(‘ [0 Hz] ‘,f1) +
     format(‘A2= %.2n’,[a2])+ formatfloat(‘ [0 Hz]’,f2);
 memo1.Lines.add(s);
 chart1.Title.Text[0]:= ‘Декодирован ‘ +s
end;
procedure TForm1.Edit1KeyPress(Sender: TObject; var Key: Char);
begin
 if key=#13 then dtmf1.DTMF_keys:= edit1.Text
end;
procedure TForm1.DTMF1Spektra(Sender: TObject; series: TBarSeries);
begin
  series1.Assign(series)
end;Syhi-подсветка кода

После чего, запустив проект на компиляцию, введем строку для генерации. Нажав ENTER можем наблюдать спектр сигнала и декодированные посылки (см. рисунок 8):

Рис. 8. Окно тестового модуля спектроанализатора-кодера-декодера DTMF

Данная методика, с добавлением анализатора файлов WAV/MP3, была использована при разработке многофункционального ПО «SPEKTRA» [7, 8] (см. рисунок 9). Следует обратить ваше внимание, что разработанный декодер DTMF не полностью соответствует всем требованиям ITU-T, т.к. не учтен анализ вторичных гармоник, например при речевом сигнале. Что это означает? Как вы наверняка знаете, человеческая речь и речь вообще, музыка характеризуются большим количеством гармоник, а сам DTMF сигнал обладает низким уровнем его вторичных гармоник. Поэтому для избежания ложного детектирования DTMF сигнала, необходимо добавить оценку этого уровня и сравнивать его с уровнями определенных частот в спектре для речи и музыки. Эти характерные частоты (форманты)** для речи в принципе известны, например характерный мужской бас сосредоточен в области 200-250 Гц и женский голос ближе к 500-700 Гц.

Но это уже тема для отдельной статьи…

Рис. 9. Отсчеты и спектр сигнала DTMF. Буква «A»

Заключение

Полные исходные тексты компонента спектроанализатора-кодера-декодера DTMF и ресурсы тестового проекта приведены в виде ресурсов в теме «Журнал клуба программистов. Восьмой выпуск» или непосредственно в архиве с журналом.

Ресурсы

AVR314: Двутональный DTMF генератор http://www.gaw.ru/html.cgi/txt/app/micros/avr/index.htm
MSP430: DTMF-Controlled http://www.gaw.ru/html.cgi/txt/app/micros/msp430/index.htm
Practical Design Techniques for Sensor Signal Conditioning, Analog Devices, 1998
DTMF Detector Data Sheet, 2001 http://www.miketdspsolutions.com/dtmf.pdf
Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов. – М., Мир, 1978
Л.А.Осипов. Обработка сигналов на цифровых процессорах. Линейно-аппроксимирующий метод. – М., Горячая линия — Телеком, 2001
Е.Бадло, С.Бадло. Виртуальные приборы. Спектроанализатор своими руками. – Радиолюбитель, Минск, 2009, №3, с.32-36 http://raxp.radioliga.com/cnt/s.php?p=v3.pdf или
Ресурсы анализатора спектра файлов и сигналов DTMF-FFT http://raxp.radioliga.com/cnt/s.php?p=dtmf_res.zip
Ресурсы компонента спектроанализатора-кодера-декодера DTMF http://raxp.radioliga.com/cnt/s.php?p=fft.zip

Статья из восьмого выпуска журнала «ПРОграммист».

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