Блог программистов

Автор: Бадло Сергей Григорьевич aka raxp
http://raxp.radioliga.com

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

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

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

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

* В 1822 году Фурье опубликовал работу «Аналитическая теория тепла», сыгравшую значительную роль в последующей истории математики. В этой работе он описал метод разделения переменных (преобразование Фурье), основанный на представлении функций тригонометрическими рядами (ряды Фурье). Фурье также сделал попытку доказать возможность разложения в тригонометрический ряд любой произвольной функции и, хоть его попытка оказалась неудачна, она фактически стала основой современной цифровой обработки сигналов. Как известно, звуковой сигнал является суперпозицией звуковых колебаний различных частот, то есть такой сигнал можно представить в виде «классического» спектра, описываемого f(?). Именно преобразование Фурье однозначно определяет соответствие между j(t) и f(?) и лежит в основе Фурье-спектроскопии / Лит.

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

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

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

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

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

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

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

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

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

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

Период колебания равен 2pi. Если мы хотим увеличить период до T, то надо умножить переменную t на коэффициент. Это вызовет растяжение графика по горизонтали: 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;

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

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

Разрабатываемый компонент предназначен для построения спектра аудио-сигнала, кодирования и декодирования двух-тоновых посылок 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

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

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

• среда разработки Borland Delphi 5-7
• аудиокарта

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

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

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

ЛИСТИНГ-1
основной модуль
…
unit DTMF;interfaceuses MMSystem, Windows, SysUtils, Messages, Classes, controls, extctrls, series, TeEngine, math;

type // тип данных wave- ind
TData16           = array [0..127] of smallint;
PData16           = ^TData16;

Type // установки для waveform
SINEWAVE = packed record
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) of object; // выдача DTMF
Tspektr = procedure(Sender:TObject; series: TbarSeries) of object; // выдача спектра
TdataOsc = procedure(Sender:TObject; series: TfastlineSeries) of object; // сырой набор данных

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                  : string read Fabout write setabout;
property DTMF_keys           : string read Fkey write gen_dtmf;
property DTMF_volume        : integer read Fvol write fvol default 100;
property DTMF_duration_ms : integer read Flen write flen default 250;
property FFT_point             : integer read 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;

procedure Register;

const // таблица соответствия частот DTMF
keys   = '1234567890*#abcd';
DTMF1: array [1..16] of integer
=(697,697,697,770,770,770,852,852,852,941,941,941,697,770,852,941);
DTMF2: array [1..16] of integer
=(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;

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;

implementation

Обычно, генерацию звука в памяти и воспроизведение в среде 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 do begin
dblAmplitudeL:= 0.25;
dblFrequency := freq; // задаем частоты генерации WAVEFORM
dblVolumeF := dblVolume
end
end;

intSineCount := length(SineWaves)-1;
for i:=0 to intSineCount do begin
dblWaveTime := 1 / SineWaves.dblFrequency;
dblSampleTime := 1 / lngSampleRate;
SineWaves.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 := 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= lowercase(fkey[gl_key]) then begin
Freq[0]:= dtmf1; // задаем частоты-
Freq[1]:= dtmf2;
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 -----------------------------------------------------

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


//------------------------------------
// ДЕКОДЕР DTMF + СПЕКТРОАНАЛИЗ
//------------------------------------
function FFT(var x, y:array of double;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 div 2;
while rt<nn do
rt:=rt*2;nn:=rt*2;
z:=2*pi/nn;case nf of  // выбор окна подавления (Блэкмана-Хэрриса)-
//---~67dB-----
1:begin
a0:=0.42323;
a1:=0.49755;
a2:=0.07922;
a3:=0;
end;
2:begin
a0:=0.44959;
a1:=0.49364;
a2:=0.05677;
a3:=0;
end;
//---~92dB-----
3:begin
a0:=0.35875;
a1:=0.48829;
a2:=0.14128;
a3:=0.01168;
end;
//---~72dB-----
4:begin
a0:=0.40217;
a1:=0.49703;
a2:=0.09392;
a3:=0.00183;
end;
5:begin // без изменений-
a0:=1;
a1:=0;
a2:=0;
a3:=0;
end
end;

for i:=0  to nn-1 do begin
w:=a0-a1*cos(z*i)+a2*cos(z*2*i)+a3*cos(z*3*i);
x:=x*w;
y:=y*w;
end;
//-------------------------------------------

ll:= nn;
M := nn div 2;
M1:= Nn-1;
while ll>=2 do begin
l:=ll div 2;
u1:=1;
u2:=0;
t1:=PI/l;
c:=cos(t1);
s:=(-1)*ii*sin(t1);
for j:=0  to l-1 do
begin
i:=j;
while i<nn do
begin
p:=i+l;
t1:=x+x[p];
t2:=y+y[p];
t3:=x-x[p];
t4:=y-y[p];
x[p]:=t3*u1-t4*u2;
y[p]:=t4*u1+t3*u2;
x:=t1;
y:=t2;
i:=i+ll
end;

u3:=u1*c-u2*s;
u2:=u2*c+u1*s;
u1:=u3
end;

ll:=ll div 2
end;

j:=0;
for i:=0  to m1-1 do begin
if i>j then begin
t1:=x[j];
t2:=y[j];
x[j]:=x;
y[j]:=y;
x:=t1;
y:=t2
end;

k:=m;
while j>=k do begin
j:=j-k;
k:=k div 2
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:=  Seriesin.YValue;
b:= 0
end;

FFT(a, b, k, f_window, 1);                   // домножаем на выбранное окно, получение спектра
for i:=0  to k div 2-1 do begin               // отсекаем зеркалку
d:= sqrt(a*a + b*b);              // получаем модуль из 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^)
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;
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;
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*0.98<f1) and (DTMF2*1.02>f1) and  // 1 амплитуда >2
(DTMF1*0.98<f2) and (DTMF1*1.02>f2) then begin
signal:= keys;
break
end;

if (DTMF1*0.98<f1) and (DTMF1*1.02>f1) and // 1 амплитуда >2
(DTMF2*0.98<f2) and (DTMF2*1.02>f2) then begin
signal:= keys;
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; // без сглаживания
dB_54: f_window:= 2; // по уровню 54 дБ
dB_67: f_window:= 1; // 67 дБ
dB_72: f_window:= 4; // 72 дБ
dB_92: f_window:= 3  // 92 дБ
end
end;

procedure TDTMF.setabout(const Value: string);
begin
fabout:= 'Badlo Sergey'
end;

//============================================================
procedure Register;
begin
RegisterComponents('RAMEDIA', [TDTMF])
end;

end.
…

Осталось проверить работоспособность нашего модуля. Для этого, выбрав в меню «Component/Install Component» проинсталлируем компонент DTMF (модуль <dtmf.pas>) (см. рис.5-7):

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

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

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

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

ЛИСТИНГ-2
тестовый модуль
…
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;
…

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

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

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

** именно за счет использования этой особенности речи в телефонной и радиосвязи полоса пропускания канала ограничивается от 300 до 3400 Гц, что позволяет применить различные компандеры, системы частотного мультиплексирования. Причем, эти частоты приняты МСЭ-Т в качестве границ эффективного спектра речи.

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

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

Заключение

Полные исходные тексты компонента спектроанализатора-кодера-декодера DTMF и ресурсы тестового проекта (файл fft.zip) вы можете загрузить тут на форуме или с сайта автора http://raxp.radioliga.com [9]

Если тема представляет для вас интерес — пишите, задавайте вопросы на форуме клуба программистов http://www.programmersforum.ru

Ресурсы

• 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 или http://www.programmersforum.ru/attachment.php?attachmentid=20723&d=1264659713
• Ресурсы анализатора спектра файлов и сигналов DTMF-FFT http://raxp.radioliga.com/cnt/s.php?p=dtmf_res.zip
• Ресурсы компонента спектроанализатора-кодера-декодера DTMF http://raxp.radioliga.com/cnt/s.php?p=fft.zip