Воскресенье, 28 Апр 2024, 14:00
Uchi.ucoz.ru
Меню сайта
Форма входа

Категории раздела
Высшая математика [11]
Экономическая социология [95]
Основы Менеджмента [64]
Бухгалтерский учёт [157]
Философия [163]
Мировая Экономика [603]
Бизнес планирование [29]
Финансирование и кредитование инвест [105]
Ценообразование [46]
Гражданское право [196]
Права Человека [173]
Основы Маркетинга [207]
Основы энергосбережения [55]
Информатика [0]
Экология и устойчивое развитие [0]
Физика для студентов [0]
Основы права [0]
Политология [0]
Не стандартные примеры на Delphi [169]
Примеры на Delphi7 [108]
Алгоритмы [94]
API [110]
Pascal [152]
Базы Данных [6]
Новости
Чего не хватает сайту?
500
Статистика
Зарегистрировано на сайте:
Всего: 51635

Онлайн всего: 2
Гостей: 2
Пользователей: 0
Яндекс.Метрика
Рейтинг@Mail.ru

Каталог статей

Главная » Статьи » Студентам » Алгоритмы
Действительно БЫСТРОЕ преобразование сигнала в спекр и обратно (методы Хартли, Фурье и классический)
Публикую присланный читателем алгоритм:

{$A+,B-,C+,D+,E-,F-,G+,H+,I+,J+,K-,L+,M-,N+,O-,P+,Q-,R-,S-,T-,U-,V+,W-,X+,Y+
,Z1}

{$MINSTACKSIZE
$00004000}

{$MAXSTACKSIZE $00100000}

{$IMAGEBASE
$00400000}

{$APPTYPE GUI}

unit
Main;

interface

uses

Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms,
Dialogs,
Buttons, ExtCtrls, ComCtrls, Menus;
type

TfmMain = class(TForm)
MainMenu1: TMainMenu;
N1: TMenuItem;
N2: TMenuItem;
StatusBar1: TStatusBar;
N3: TMenuItem;
imgInfo: TImage;
Panel1: TPanel;
btnStart: TSpeedButton;
procedure btnStartClick(Sender: TObject);
procedure FormCreate(Sender: TObject);
procedure FormClose(Sender: TObject; var Action:
TCloseAction);
end;
var

fmMain: TfmMain;
implementation

Uses
PFiles;

{$R *.DFM}

function
Power2(lPower: Byte): LongInt;

begin
Result := 1 Shl lPower;
end;

procedure ClassicDirect(Var
aSignal, aSpR, aSpI: Array Of Double; N:
LongInt);

var lSrch : LongInt;
var lGarm : LongInt;
var dSumR : Double;
var dSumI : Double;
begin
for lGarm := 0 to N div
2 - 1 do begin
dSumR := 0;
dSumI := 0;
for lSrch := 0 to N - 1 do begin
dSumR := dSumR + aSignal[lSrch] * Cos(lGarm * lSrch / N * 2 * PI);
dSumI := dSumI + aSignal[lSrch] * Sin(lGarm * lSrch / N * 2 * PI);
end;
aSpR[lGarm] := dSumR;
aSpI[lGarm] := dSumI;
end;
end;

procedure
ClassicInverce(Var aSpR, aSpI, aSignal: Array Of Double;
N:
LongInt);

var lSrch : LongInt;
var lGarm : LongInt;
var dSum : Double;
begin
for lSrch := 0 to N-1 do begin
dSum := 0;
for lGarm := 0 to N
div 2 -1
do dSum := dSum
+ aSpR[lGarm] * Cos(lSrch * lGarm * 2 * Pi
/ N)
+ aSpI[lGarm] * Sin(lSrch * lGarm * 2 * Pi
/ N);
aSignal[lSrch] := dSum*2;
end;
end;

Function InvertBits(BF, DataSize,
Power: Word) : Word;

Var BR : Word;
Var NN : Word;
Var L : Word;
Begin
br:= 0;
nn:= DataSize;
For l:= 1 To Power
Do
Begin
NN:= NN Div 2;
If (BF >= NN) Then
Begin
BR:= BR + Power2(l - 1);
BF:= BF - NN
End;
End;
InvertBits:=BR;
End;
Procedure FourierDirect(Var
RealData,VirtData,ResultR,ResultV: Array Of
Double; DataSize:
LongInt);

Var A1 : Real;
Var A2 : Real;
Var B1 : Real;
Var B2 : Real;
Var D2 : Word;
Var C2 : Word;
Var C1 : Word;
Var D1 : Word;
Var I : Word;
Var J : Word;
Var K : Word;
Var Cosin : Real;
Var Sinus : Real;
Var wIndex : Word;
Var Power : Word;
Begin
C1:= DataSize Shr 1;
C2:= 1;
for Power:=0 to 15 //hope it will be faster
thenround(ln(DataSize)/ln(2))

do if Power2(Power)=DataSize
then Break;
For I:= 1 To Power Do
Begin
D1:= 0;
D2:= C1;
For J:= 1 To C2 Do
Begin
wIndex:=InvertBits(D1 Div C1, DataSize, Power);
Cosin:= +(Cos((2 * Pi /
DataSize)*wIndex));
Sinus:= -(Sin((2 * Pi /
DataSize)*wIndex));
For K:= D1 To D2 - 1 Do
Begin
A1:= RealData[K];
A2:= VirtData[K];
B1:= ((Cosin * RealData[K + C1] - Sinus * VirtData[K + C1])
);
B2:= ((Sinus * RealData[K + C1] + Cosin * VirtData[K + C1])
);
RealData[K]:= A1 + B1 ;
VirtData[K]:= A2 + B2 ;
RealData[K + C1]:= A1 - B1;
VirtData[K + C1]:= A2 - B2;
End;
Inc(D1,C1 * 2);
Inc(D2,C1 * 2);
End;
C1:=C1 Div 2;
C2:=C2 * 2;
End;
For I:= 0 To DataSize
Div 2 -1 Do
Begin
ResultR[I]:= + RealData[InvertBits(I, DataSize, Power)];
ResultV[I]:= - VirtData[InvertBits(I, DataSize, Power)];
End;
End;

Procedure Hartley(iSize:
LongInt;Var aData : Array Of Double);

Type taDouble = Array[0..MaxLongInt Div SizeOf(Double)-1] Of Double;
Var prFI,prFN,prGI : ^taDouble;
Var rCos,rSin : Double;
Var rA,rB,rTemp : Double;
Var rC1,rC2,rC3,rC4 : Double;
Var rS1,rS2,rS3,rS4 : Double;
Var rF0,rF1,rF2,rF3 : Double;
Var rG0,rG1,rG2,rG3 : Double;
Var iK1,iK2,iK3,iK4 : LongInt;
Var iSrch,iK,iKX : LongInt;
Begin
iK2:=0;
For iK1:=1 To iSize-1 Do Begin
iK:=iSize Shr 1;
Repeat
iK2:=iK2 Xor iK;
If (iK2 And iK)<>0
Then Break;
iK:=iK Shr 1;
Until False;
If iK1>iK2 Then Begin
rTemp:=aData[iK1];
aData[iK1]:=aData[iK2];
aData[iK2]:=rTemp;
End;
End;
iK:=0;
While (1 Shl iK)<iSize
Do Inc(iK);
iK:=iK And 1;
If iK=0 Then Begin
prFI:=@aData;
prFN:=@aData;
prFN := @prFN[iSize];
While Word(prFI)<Word(prFN) Do Begin
rF1:=prFI^[0]-prFI^[1];
rF0:=prFI^[0]+prFI^[1];
rF3:=prFI^[2]-prFI^[3];
rF2:=prFI^[2]+prFI^[3];
prFI^[2]:=rF0-rF2;
prFI^[0]:=rF0+rF2;
prFI^[3]:=rF1-rF3;
prFI^[1]:=rF1+rF3;
prFI := @prFI[4];
End;
End Else Begin
prFI:=@aData;
prFN:=@aData;
prFN := @prFN[iSize];
prGI:=prFI;
prGI := @prGI[1];
While Word(prFI)<Word(prFN) Do begin
rC1:=prFI^[0]-prGI^[0];
rS1:=prFI^[0]+prGI^[0];
rC2:=prFI^[2]-prGI^[2];
rS2:=prFI^[2]+prGI^[2];
rC3:=prFI^[4]-prGI^[4];
rS3:=prFI^[4]+prGI^[4];
rC4:=prFI^[6]-prGI^[6];
rS4:=prFI^[6]+prGI^[6];
rF1:=rS1-rS2;
rF0:=rS1+rS2;
rG1:=rC1-rC2;
rG0:=rC1+rC2;
rF3:=rS3-rS4;
rF2:=rS3+rS4;
rG3:=Sqrt(2)*rC4;
rG2:=Sqrt(2)*rC3;
prFI^[4]:=rF0-rF2;
prFI^[0]:=rF0+rF2;
prFI^[6]:=rF1-rF3;
prFI^[2]:=rF1+rF3;
prGI^[4]:=rG0-rG2;
prGI^[0]:=rG0+rG2;
prGI^[6]:=rG1-rG3;
prGI^[2]:=rG1+rG3;
prFI := @prFI[8];
prGI := @prGI[8];
End;
End;
If iSize<16 Then
Exit;
Repeat
Inc(iK,2);
iK1:=1 Shl iK;
iK2:=iK1 Shl 1;
iK4:=iK2 Shl 1;
iK3:=iK2+iK1;
iKX:=iK1 Shr 1;
prFI:=@aData;
prGI:=prFI;
prGI := @prGI[iKX];
prFN:=@aData;
prFN := @prFN[iSize];
Repeat
rF1:= prFI^[000]-prFI^[iK1];
rF0:= prFI^[000]+prFI^[iK1];
rF3:= prFI^[iK2]-prFI^[iK3];
rF2:= prFI^[iK2]+prFI^[iK3];
prFI^[iK2]:=rF0-rF2;
prFI^[000]:=rF0+rF2;
prFI^[iK3]:=rF1-rF3;
prFI^[iK1]:=rF1+rF3;
rG1:=prGI^[0]-prGI^[iK1];
rG0:=prGI^[0]+prGI^[iK1];
rG3:=Sqrt(2)*prGI^[iK3];
rG2:=Sqrt(2)*prGI^[iK2];
prGI^[iK2]:=rG0-rG2;
prGI^[000]:=rG0+rG2;
prGI^[iK3]:=rG1-rG3;
prGI^[iK1]:=rG1+rG3;
prGI := @prGI[iK4];
prFI := @prFI[iK4];
Until Not (Word(prFI)<Word(prFN));
rCos:=Cos(Pi/2/Power2(iK));
rSin:=Sin(Pi/2/Power2(iK));
rC1:=1;
rS1:=0;
For iSrch:=1 To iKX-1 Do Begin
rTemp:=rC1;
rC1:=(rTemp*rCos-rS1*rSin);
rS1:=(rTemp*rSin+rS1*rCos);
rC2:=(rC1*rC1-rS1*rS1);
rS2:=(2*(rC1*rS1));
prFN:=@aData;
prFN := @prFN[iSize];
prFI:=@aData;
prFI := @prFI[iSrch];
prGI:=@aData;
prGI := @prGI[iK1-iSrch];
Repeat
rB:=(rS2*prFI^[iK1]-rC2*prGI^[iK1]);
rA:=(rC2*prFI^[iK1]+rS2*prGI^[iK1]);
rF1:=prFI^[0]-rA;
rF0:=prFI^[0]+rA;
rG1:=prGI^[0]-rB;
rG0:=prGI^[0]+rB;
rB:=(rS2*prFI^[iK3]-rC2*prGI^[iK3]);
rA:=(rC2*prFI^[iK3]+rS2*prGI^[iK3]);
rF3:=prFI^[iK2]-rA;
rF2:=prFI^[iK2]+rA;
rG3:=prGI^[iK2]-rB;
rG2:=prGI^[iK2]+rB;
rB:=(rS1*rF2-rC1*rG3);
rA:=(rC1*rF2+rS1*rG3);
prFI^[iK2]:=rF0-rA;
prFI^[0]:=rF0+rA;
prGI^[iK3]:=rG1-rB;
prGI^[iK1]:=rG1+rB;
rB:=(rC1*rG2-rS1*rF3);
rA:=(rS1*rG2+rC1*rF3);
prGI^[iK2]:=rG0-rA;
prGI^[0]:=rG0+rA;
prFI^[iK3]:=rF1-rB;
prFI^[iK1]:=rF1+rB;
prGI := @prGI[iK4];
prFI := @prFI[iK4];
Until Not (LongInt(prFI) < LongInt(prFN));
End;
Until Not (iK4<iSize);
End;

Procedure
HartleyDirect(
Var aData : Array Of
Double;

iSize : LongInt);Var rA,rB
: Double;
Var iI,iJ,iK : LongInt;

Begin
Hartley(iSize,aData);
iJ:=iSize-1;
iK:=iSize Div 2;
For iI:=1 To iK-1 Do Begin
rA:=aData[ii];
rB:=aData[ij];
aData[iJ]:=(rA-rB)/2;
aData[iI]:=(rA+rB)/2;
Dec(iJ);
End;
End;
Procedure HartleyInverce(
Var
aData : Array Of Double;

iSize : LongInt);

Var rA,rB : Double;
Var iI,iJ,iK : LongInt;
Begin
iJ:=iSize-1;
iK:=iSize Div 2;
For iI:=1 To iK-1 Do Begin
rA:=aData[iI];
rB:=aData[iJ];
aData[iJ]:=rA-rB;
aData[iI]:=rA+rB;
Dec(iJ);
End;
Hartley(iSize,aData);
End;
//not
tested
procedure HartleyDirectComplex(real,imag: Array
of Double;n: LongInt);
var a,b,c,d :
double;

q,r,s,t : double;
i,j,k : LongInt;begin

j:=n-1;
k:=n div 2;
for i:=1 to k-1 do begin
a := real[i]; b := real[j]; q:=a+b; r:=a-b;
c := imag[i]; d := imag[j]; s:=c+d; t:=c-d;
real[i] := (q+t)*0.5; real[j] :=
(q-t)*0.5;
imag[i] := (s-r)*0.5; imag[j] :=
(s+r)*0.5;
dec(j);
end;
Hartley(N,Real);
Hartley(N,Imag);end;

//not tested
procedure
HartleyInverceComplex(real,imag: Array Of Double;N:
LongInt);
var a,b,c,d :double;

q,r,s,t :double;
i,j,k :longInt;
begin
Hartley(N,real);
Hartley(N,imag);
j:=n-1;
k:=n div 2;
for i:=1 to k-1 do begin
a := real[i]; b := real[j]; q:=a+b; r:=a-b;
c := imag[i]; d := imag[j]; s:=c+d; t:=c-d;
imag[i] := (s+r)*0.5; imag[j] :=
(s-r)*0.5;
real[i] := (q-t)*0.5; real[j] :=
(q+t)*0.5;
dec(j);
end;
end;

procedure DrawSignal(var
aSignal: Array Of Double;N,lColor : LongInt);

var lSrch : LongInt;
var lHalfHeight : LongInt;
begin
with fmMain do begin
lHalfHeight := imgInfo.Height Div 2;
imgInfo.Canvas.MoveTo(0,lHalfHeight);
imgInfo.Canvas.Pen.Color := lColor;
for lSrch := 0 to N-1 do begin
imgInfo.Canvas.LineTo(lSrch,Round(aSignal[lSrch]) +
lHalfHeight);
end;
imgInfo.Repaint;
end;
end;

procedure DrawSpector(var
aSpR, aSpI: Array Of Double;N, lColR, lColI :
LongInt);

var lSrch : LongInt;
var lHalfHeight : LongInt;
begin
with fmMain do begin
lHalfHeight := imgInfo.Height Div 2;
for lSrch := 0 to N
Div 2 do begin
imgInfo.Canvas.Pixels[lSrch ,Round(aSpR[lSrch]/N) + lHalfHeight]
:=lColR;

imgInfo.Canvas.Pixels[lSrch + N Div 2 ,Round(aSpI[lSrch]/N) +lHalfHeight] := lColI;

end;
imgInfo.Repaint;
end;
end;
const N =
512;
var aSignalR :
Array[0..N-1] Of
Double;//
var aSignalI :
Array[0..N-1] Of
Double;//
var aSpR, aSpI :
Array[0..N Div 2-1] Of Double;//
var lFH :
LongInt;

procedure TfmMain.btnStartClick(Sender: TObject);

const Epsilon = 0.00001;
var lSrch : LongInt;
var aBuff : Array[0..N-1] Of ShortInt;
begin
if lFH > 0 then
begin// Repeat

if F.Read(lFH,@aBuff,N) <> N then begin
Exit;
end;
for lSrch := 0 to N-1 do begin
aSignalR[lSrch] := ShortInt(aBuff[lSrch]+$80);
aSignalI[lSrch] := 0;
end;

imgInfo.Canvas.Rectangle(0,0,imgInfo.Width,imgInfo.Height);
DrawSignal(aSignalR, N, $D0D0D0);
//
ClassicDirect(aSignalR, aSpR, aSpI, N); //result in aSpR &
aSpI,
aSignal unchanged
// FourierDirect(aSignalR, aSignalI, aSpR,
aSpI, N); //result in aSpR &
aSpI, aSiggnalR & aSignalI
modified

HartleyDirect(aSignalR, N); //result in source aSignal ;-)

DrawSpector(aSignalR, aSignalR[N Div 2 -1], N, $80, $8000);
DrawSpector(aSpR, aSpI, N, $80, $8000);
{ for lSrch := 0 to N
div 2 -1 do begin //comparing classic & Hartley

if (Abs(aSpR[lSrch] - aSignal[lSrch]) > Epsilon)
or ((lSrch > 0) And (Abs(aSpI[lSrch] - aSignal[N - lSrch]) >
Epsilon))
then MessageDlg('Error comparing',mtError,[mbOK],-1);
end;}

HartleyInverce(aSignalR, N); //to restore original signal withHartleyDirect
//
ClassicInverce(aSpR, aSpI, aSignalR, N); //to restore
original
signal with ClassicDirect or FourierDirect

for lSrch := 0 to N -1
do aSignalR[lSrch]:= aSignalR[lSrch]/N;
//scaling

DrawSignal(aSignalR, N, $D00000);
Application.ProcessMessages;// Until
False;

end;
end;
procedure TfmMain.FormCreate(Sender:
TObject);

begin
lFH := F.Open('input.pcm',
ForRead);
end;
procedure TfmMain.FormClose(Sender:
TObject; var Action: TCloseAction);

begin
F.Close(lFH);
end;
end.
Категория: Алгоритмы | Добавил: Lerka (21 Ноя 2012)
Просмотров: 475 | Рейтинг: 1.0/ 4 Оштрафовать | Жаловаться на материал
Похожие материалы
Всего комментариев: 0

Для блога (HTML)


Для форума (BB-Code)


Прямая ссылка

Профиль
Воскресенье
28 Апр 2024
14:00

Вы из группы: Гости
Вы уже дней на сайте
У вас: непрочитанных сообщений
Добавить статью
Прочитать сообщения
Регистрация
Вход
Улучшенный поиск
Поиск по сайту Поиск по всему интернету
Наши партнеры
Интересное
Популярное статьи
Портфолио ученика начальной школы
УХОД ЗА ВОЛОСАМИ ОЧЕНЬ ПРОСТ — ХОЧУ Я ЭТИМ ПОДЕЛИТ...
Диктанты 2 класс
Детство Л.Н. Толстого
Библиографический обзор литературы о музыке
Авторская программа элективного курса "Практи...
Контрольная работа по теме «Углеводороды»
Поиск
Главная страница
Используются технологии uCoz