1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258
| // Gunnar Bolle, FPrefect@t-online.de
// Simple FFT component
// Feel free to use this code
// Based upon a pascal routine from - Don Cross <dcross@intersrv.com>
// Thanks Don, nice Job.
// If you want to know more about FFT and its theorie visit
// http://www.intersrv.com/~dcross
// Thanks to Kees Huls for providing me the missing author information
//
// Please note : I didn't have the time to write a sample application for this.
// Please do not ask me about how to use this one ...
// If you're aware of FFT you'll know how. Otherwise, try
// some of those neat Button components at DSP. They're quite easy
// to handle.
unit DSXFastFourier;
interface
uses
Windows, Messages, SysUtils, Classes,math;
procedure Register;
type
TComplex = Record
Real : double;
imag : double;
end;
TOnGetDataEvent = procedure(index : integer; var Value : TComplex) of Object;
TComplexArray = array [0..0] of TComplex;
PComplexArray = ^TComplexArray;
EFastFourierError = class(Exception);
TDSXFastFourier = Class(TComponent)
private
FNumSamples : integer;
FInBuffer : PComplexArray;
FOutBuffer : PComplexArray;
FOnGetData : TOnGetDataEvent;
function IsPowerOfTwo ( x: word ): boolean;
function NumberOfBitsNeeded ( PowerOfTwo: word ): word;
function ReverseBits ( index, NumBits: word ): word;
procedure FourierTransform ( AngleNumerator: double );
procedure SetNumSamples(value : integer);
function GetTransformedData(idx : integer) : TComplex;
constructor
create(AOwner : TComponent);
destructor
destroy;
public
procedure fft;
procedure ifft;
procedure CalcFrequency (FrequencyIndex: word);
published
property OnGetData : TOnGetDataEvent read FOnGetData write FOnGetData;
property NumSamples : integer read FNumSamples write SetNumSamples;
property SampleCount : Integer read FNumSamples;
property TransformedData[idx : integer] : TComplex read GetTransformedData;
end;
implementation
constructor TDSXFastFourier.Create(AOwner : TComponent);
begin
inherited create(AOwner);
end;
destructor TDSXFastFourier.Destroy;
begin
if Assigned(FInBuffer) then
FreeMem(FinBuffer);
if Assigned(FOutBuffer) then
FreeMem(FOutBuffer);
end;
procedure TDSXFastFourier.SetNumSamples(value : integer);
begin
FNumSamples := value;
if Assigned(FInBuffer) then
FreeMem(FinBuffer);
if Assigned(FOutBuffer) then
FreeMem(FOutBuffer);
try
getMem(FInBuffer, sizeof(TComplex)*FNumSamples);
getMem(FOutBuffer, sizeof(TComplex)*FNumSamples);
except on EOutOfMemory do
raise EFastFourierError.Create('Could not allocate memory for complex arrays');
end;
end;
function TDSXFastFourier.GetTransformedData(idx : integer) : TComplex;
begin
Result := FOutBuffer[idx];
end;
function TDSXFastFourier.IsPowerOfTwo ( x: word ): boolean;
var i, y: word;
begin
y := 2;
for i := 1 to 31 do begin
if x = y then begin
IsPowerOfTwo := TRUE;
exit;
end;
y := y SHL 1;
end;
IsPowerOfTwo := FALSE;
end;
function TDSXFastFourier.NumberOfBitsNeeded ( PowerOfTwo: word ): word;
var i: word;
begin
for i := 0 to 16 do begin
if (PowerOfTwo AND (1 SHL i)) <> 0 then begin
NumberOfBitsNeeded := i;
exit;
end;
end;
end;
function TDSXFastFourier.ReverseBits ( index, NumBits: word ): word;
var i, rev: word;
begin
rev := 0;
for i := 0 to NumBits-1 do begin
rev := (rev SHL 1) OR (index AND 1);
index := index SHR 1;
end;
ReverseBits := rev;
end;
procedure TDSXFastFourier.FourierTransform ( AngleNumerator: double);
var
NumBits, i, j, k, n, BlockSize, BlockEnd: word;
delta_angle, delta_ar: double;
alpha, beta: double;
tr, ti, ar, ai: double;
begin
if not IsPowerOfTwo(FNumSamples) or (FNumSamples<2) then
raise EFastFourierError.Create('NumSamples is not a positive integer power of 2');
if not assigned(FOnGetData) then
raise EFastFourierError.Create('You must specify an OnGetData handler');
NumBits := NumberOfBitsNeeded (FNumSamples);
for i := 0 to FNumSamples-1 do begin
j := ReverseBits ( i, NumBits );
FOnGetData(i,FInBuffer[i]);
FOutBuffer[j] := FInBuffer[i];
end;
BlockEnd := 1;
BlockSize := 2;
while BlockSize <= FNumSamples do begin
delta_angle := AngleNumerator / BlockSize;
alpha := sin ( 0.5 * delta_angle );
alpha := 2.0 * alpha * alpha;
beta := sin ( delta_angle );
i := 0;
while i < FNumSamples do begin
ar := 1.0; (* cos(0) *)
ai := 0.0; (* sin(0) *)
j := i;
for n := 0 to BlockEnd-1 do begin
k := j + BlockEnd;
tr := ar*FOutBuffer[k].Real - ai*FOutBuffer[k].Imag;
ti := ar*FOutBuffer[k].Imag + ai*FOutBuffer[k].Real;
FOutBuffer[k].Real := FOutBuffer[j].Real - tr;
FOutBuffer[k].Imag := FOutBuffer[j].Imag - ti;
FOutBuffer[j].Real := FOutBuffer[j].Real + tr;
FOutBuffer[j].Imag := FOutBuffer[j].Imag + ti;
delta_ar := alpha*ar + beta*ai;
ai := ai - (alpha*ai - beta*ar);
ar := ar - delta_ar;
INC(j);
end;
i := i + BlockSize;
end;
BlockEnd := BlockSize;
BlockSize := BlockSize SHL 1;
end;
end;
procedure TDSXFastFourier.fft;
begin
FourierTransform ( 2*PI);
end;
procedure TDSXFastFourier.ifft;
var
i: word;
begin
FourierTransform ( -2*PI);
(* Normalize the resulting time samples... *)
for i := 0 to FNumSamples-1 do begin
FOutBuffer[i].Real := FOutBuffer[i].Real / FNumSamples;
FOutBuffer[i].Imag := FOutBuffer[i].Imag / FNumSamples;
end;
end;
procedure TDSXFastFourier.CalcFrequency (FrequencyIndex: word);
var
k: word;
cos1, cos2, cos3, theta, beta: double;
sin1, sin2, sin3: double;
begin
FOutBuffer[0].Real := 0.0;
FOutBuffer[0].Imag := 0.0;
theta := 2*PI * FrequencyIndex / FNumSamples;
sin1 := sin ( -2 * theta );
sin2 := sin ( -theta );
cos1 := cos ( -2 * theta );
cos2 := cos ( -theta );
beta := 2 * cos2;
for k := 0 to FNumSamples-1 do begin
sin3 := beta*sin2 - sin1;
sin1 := sin2;
sin2 := sin3;
cos3 := beta*cos2 - cos1;
cos1 := cos2;
cos2 := cos3;
FOutBuffer[0].Real := FOutBuffer[0].Real + FInBuffer[k].Real*cos3 - FInBuffer[k].Imag*sin3;
FOutBuffer[0].Imag := FOutBuffer[0].Imag + FInBuffer[k].Imag*cos3 + FInBuffer[k].Real*sin3;
end;
end;
procedure Register;
begin
RegisterComponents('Free', [TDSXFastFourier]);
end;
end. |
Partager