» Быстрое преобразование Фурье. Практика использования Borland Delphi. Новости. . Блог программистов


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




20105 Feb

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

Автор: Бадло Сергей Григорьевич 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(N2), требуемых для прямого вычисления ДПФ.

Для чего нужно быстрое преобразование Фурье? Допустим у нас есть периодическая функция изменяющаяся по закону синуса 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

Ресурсы

Комментарии

  1. Beermonza
    February 6th, 2010 | 19:31

    Материал полезный.
    Картинка в начале статьи не корректно ставит вопрос, …тут просто однозначно вот так: осциллограф ~ стофт, поскольку задета ничтожная часть частотного диапазона, не нужно провокаций, а то с дуру можно и за равенство посчитать.
    Без претензий к статье.

  2. February 6th, 2010 | 21:56

    никаких провокаций, под картинкой есть поясняющая надпись… дело в том, что этот-же метод FFT с той-же реализацией в ПО “SPEKTRA” используется и на высоких частотах, естественно аппаратно-программный http://raxp.radioliga.com/cnt/s.php?p=mv.pdf

  3. Beermonza
    February 6th, 2010 | 22:49

    Я понимаю, что где-то есть и с платой сбора данных, …но когда фигурирует WAV-формат и дискретизация 44 кГц, которой далеко не заканчивается диапазон подлежащий исследованию, то какие-то аналогии на картинке с прибором TDS3064 (600 МГц) проводить вряд ли будет корректно, даже если есть пояснения для чего ваш компонент и какие у него параметры, …ну вы понимаете, …можно и оставить, но посчитал все же обратить внимание, это именно провокация ))) …без претензий к статье и вам.

  4. February 11th, 2010 | 17:39

    Много математики, но все равно спасибо, алгоритм хороший

  5. February 15th, 2010 | 19:51

    Есть ли аналоги вашей разработке?

  6. February 18th, 2010 | 12:48

    to Android…
    математику нужно любить, а не чтобы она любила вас :)

    to Artemich…
    уточните пожалуйста ваш вопрос… тема статьи “БПФ и методы практической реализации”, распознавание DTMF и создание компонента спектроанализатора-кодера-декодера со своей стороны уже рассматриваю как обучающий и “увлекающий” материал для энтузиастов и профессионалов… аналогов, насколько мне известно, нет.

  7. March 4th, 2010 | 09:35

    А как здесь без математики? Есть другие методы изложения?

  8. Sergey
    March 21st, 2010 | 18:21

    Уважаемый raxp и все остальные помогите пожалуйста написать программу подобную этой. Сам я, к сожалению, не настолько силён в программировании. Отблагодарю. Желающие помочь пишите на майл sawd[at]yandex.ru

  9. Алексей
    March 26th, 2010 | 22:01

    спасибо за инфу, хочу написать анализатор на свой мп3 плеер. надеюсь выше написанная статья мне в этом поможет

  10. March 26th, 2010 | 23:16

    to Sergey…
    дождитесь первого выпуска журнала “Клуба программистов”, там будет статья-продолжение, как раз на вашу тему.

    to Алексей…
    - если выберете bass в качестве основы плеера, то fft там встроенный есть, не так много точек и не универсальный, но есть http://www.programmersforum.ru/showthread.php?t=89708&page=2
    - если DirectShow, то там тоже можно через DSPackDCDSPFilter… это так сказать готовые специфические технологии

  11. Дмитрий
    May 3rd, 2010 | 23:01

    У меня к вам небольшой вопрос.
    Подскажите какие переменные отвечают за Дб, а какие за Гц.
    Мне необходимо сделать чтоб я вписал допустим условие что при Х Гц и Х Дб выскакивала ошибка, такое реально сделать?

  12. May 3rd, 2010 | 23:17

    внимательней читаем комментарии в коде (см. процедуру waveInProc2()):
    =======================
    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); // отсекаем отрицательные амплитуды

    end;
    =======================

    …однако учтите, что в FFT() идет подставка в -25дБ для амплитуды (которая вам не нужна):
    =======================

    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; // приведение к дБ и нормирование

    seriesout.Add(d)
    end
    =======================

    С уважением, Сергей

  13. Дмитрий
    May 4th, 2010 | 01:12

    спасибо за помощь, уже знаю где копаться
    когда я в цикл вставил условие
    if a>30 then message(‘ХХХ’); у меня выдало ошибку((
    А можно как-то автоматизировать процесс выбора частот и амплитуд например в edite вводить? я так понимаю чтоб это сделать необходимо переменные f и а сделать глобальными и в коде самой программы сравнивать значения заданные в edite с ими?
    возможно ли с помощью этой проги узнать скважность сигнала?
    Заранее спасибо!

  14. May 4th, 2010 | 08:34

    1- глобальные, да
    2- вы должны понимать, что судить о скважности по спектру (даже не о скважности, а о ее изменении) можно только очень косвенно и то при граничных условиях, когда форма сигнала не меняется. Скважность – это отношение периода к длительности и измерять ее надо на “сырых” данных (огибающей), скажем обрабатывать тут:
    ==================================
    if Assigned(FOnDataOsc) then FOnDataOsc(self,inwav);
    ==================================
    но, как реализовать выделение и обеспечить синхронизацию – это уже ваша задача :) и не тема данной статьи.

  15. Дмитрий
    May 9th, 2010 | 23:04

    снова я, вы знаете с амплитудной составляющей у меня проблем не возникло, а вот в частотой что-то не особо, потом я вывел значения частоты и увидел что эти значения изменяются на протяжении всего диапазона 0..~40000 и поэтому я не могу подставить условие на Гц.

  16. May 10th, 2010 | 09:00

    1- не совсем понятно, как вы проверяете свое условие
    2- гармоника она никуда не девается, она стоит на месте

  17. Дмитрий
    May 17th, 2010 | 23:14

    И опять я) Вы не подскажите почему когда я эту прогу пытаюсь запустить на другом компе мне выдаёт ошибку Invalid pointer operator?

  18. May 17th, 2010 | 23:59

    …а чем отличаются компьютеры?

  19. Дмитрий
    May 18th, 2010 | 00:07

    Ну это ноутбуки уже на пяти пробовал))) незнаю что делать, тока домой прихожу на стационарнике всё запускается, там тоже запускается но на долю секунды потом ошибка((

  20. raxp
    May 18th, 2010 | 10:02

    …ну что могу сказать, под вопросом “чем отличаются компьютеры” – имелось ввиду, какое железо (аудиокарта) там и там и какие системы там и там …у себя на нетбуках и ноутбуках все ок, системы xp и семерка.

    …запустите пошаговое выполнение в среде delphi на ноутбуке и отловите где возникает ошибка :)

  21. Дмитрий
    May 18th, 2010 | 13:08

    я понял, а это может возникать из-за того что неустановленна среда программирования Delphi на других компах? я запускал скомпилированных exe-шник

  22. raxp
    May 18th, 2010 | 13:48

    нет

  23. Дмитрий
    May 18th, 2010 | 17:06

    а может проблема в антивирусах которые возможно блокируют порт?
    или наоборот у меня вирусы, хотя вроде каспер 10 стоит с последними обновлениями?

  24. Алексей
    June 9th, 2010 | 19:32

    а как получить истинное значение частоты? ведь i – это порядковый номер выборки. которая получена со звуковой карты. потом ее загоняют в fft и длина выборки сокращается в два раза. Но если сократить размер первоначальной выборки, то и диапазон охватываемых частот сократится, а это не есть ведь правильно, они же(частоты) все равно останутся, т.е. будут поступать на звуковую карту с микрофона. И как минимальную снимаемую частоту и максимальную определить?

  25. June 12th, 2010 | 13:08

    пересчет номера выборки в частоту показан выше, читайте внимательнее …насчет всего остального:
    =================
    Спектральное представление сигнала не является результатом прямых измерений. Частотные составляющие – это расчетные (с помощью БПФ) величины на основе массива выборок.

    …диапазон охватываемых частот не сокращается, осуществляется прореживание выборок :)

  26. Алексей
    June 14th, 2010 | 15:42

    Количество точек которое выводится на экран в самом конце, напрямую зависит от длины первоначальной выборки… Может я не правильно понимаю что то..
    пересчет номера выборки в частоту в этом вот участве кода?:
    d:= 20*log10(d/k + 0.000001) -25;

  27. Алексей
    June 14th, 2010 | 16:38

    Другими словами, как в строющемся диапвзоне частот узнать минимум и максимум.. даже если при сокращении выборки идет прореживание, то максимум и минимум остается неизменными по идее?!

  28. June 14th, 2010 | 20:02

    вы второй человек, не обращающий внимание на комментарии к коду :) рядом-же комментарий …мало того, уже в самих комментариях к статье вы не первый это спрашиваете:
    ====================
    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); // отсекаем отрицательные амплитуды

    end;
    ====================
    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дб это подставка, чтоб убрать фоновые шумы

    ====================

    минимум в БПФ вообще-то – это постоянная составляющая :) теперь по максимуму, СМОТРИМ комментарии в коде:

    cntval:= header.nSamplesPerSec / outwav.YValues.Count;
    for i:= 0 to (outwav.YValues.Count)-1 do begin
    a:= outwav.YValues;
    f:= i * cntval;

    очевидно, что максимальная частота определяется в процедуре инициализации захвата звукового потока – nSamplesPerSec (частотой выборки), а теперь смотрим процедуру wcard(), что мы там видим? а видим мы комментарий:
    =====================
    nSamplesPerSec:= 44100; // дискретизация, Гц
    =====================
    другой вопрос, что не следует доверять частотам выше половины частоты дискретизации и тут следует использовать накопление отсчетов …вообще, если речь об измерениях, то следует использовать внешнюю аудиокарту и никак не бытового применения, а с расширенной полосой

  29. Reqyz
    July 5th, 2010 | 09:23

    функции ExtractByte, использующаяся в вашем примере, в природе не существует, что она должна делать?

  30. Reqyz
    July 5th, 2010 | 09:46

    этот код сплош состоит из ошибок, он вообще не компилируется, возьмем хотя бы процедуру toneGenerate, как можно присваивать массиву SineWaves, не используя индексы?, и таких косяков, полный код, пожалуйста исправьте, а потом выкладывайте

  31. July 5th, 2010 | 20:40

    Уважаемый невнимальный читатель Reqyz! Обратите ваш невнимательный взор на раздел ресурсы к статье – там вы найдете полный исходный код с компиляцией, для этого и существует раздел ресурсы :)

    ===========================
    по поводу второго, вы опять недостаточно внимательны :) как вы думаете, что означает эта запись:

    type
    SINEWAVE = packed record
    dblFrequency: Double;
    dblDataSlice: Double;
    dblAmplitudeL: Double;
    dblVolumeF: Double;
    end;

    var 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];
    dblVolumeF := dblVolume[i]
    end
    end;

    …все компилируется без ошибок. Нечасто это говорю, но на этот раз скажу – учите матчасть и будет вам счастье :) Удачи.

  32. July 27th, 2010 | 06:44

    Слов нет, одни эмоции, как же я искал это преобразование недавно

  33. November 29th, 2010 | 22:35

    Хочу обратить внимание сообщества на вот этих плагиатчиков http://lfkoptima.com/home/2010-08-12-10-04-55/44-2010-08-20-06-57-16/1863-2010-10-27-12-00-08.html

    Скриншоты:
    http://img14.imageshack.us/img14/1286/screen4zf.png
    http://img522.imageshack.us/img522/5872/screen0c.png
    http://img153.imageshack.us/img153/3517/screen1ob.png
    http://img138.imageshack.us/img138/3488/screen2ym.png

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

  34. December 1st, 2010 | 17:43

    Статья прекрасна но автор Человек Мешаюший Обществу

  35. December 1st, 2010 | 17:54

    …что вы говорите, очевидно мешаю обещству тем, что бесплатно на энтузиазме занимаюсь разработками и делюсь ими с этим самым обществом :) ?

    …очевидно, общество это lfkoptima.com вот это ? даже подозреваю, что вы это и есть …неужели так трудно было сослаться на автора? Теперь, увы, на грубость ответит хостер, о чем вас (если это вы, а думаю это вы) уже предупредил.

    На сим разговор окончен. При флуде, возможность комментариев тут будет отключена. Уж с Аларом-то решим как нить данный вопрос.

  36. Dmitry
    April 24th, 2012 | 11:24

    Почему когда я запускаю, выходит окно “прекращена работа программы SPEKTRA” ?

  37. Faust
    September 4th, 2013 | 08:28

    Это у вас что-то сбилось или у меня что-то с кодировкой…. Проверти, пожалуйста формулы. Знаки вопроса, вместо символов, почему-то….

Ответить