1 | {This unit is part of United Openlibraries of Sound (uos)}
|
---|
2 | {
|
---|
3 | audacity which this file was converted from is GPLv2
|
---|
4 | http://www.audacityteam.org/about/license/
|
---|
5 | Converted By Andrew Haines
|
---|
6 | License : modified LGPL.
|
---|
7 | Fred van Stappen fiens@hotmail.com
|
---|
8 | }
|
---|
9 |
|
---|
10 | unit uos_dsp_noiseremoval;
|
---|
11 |
|
---|
12 | {$mode objfpc}{$H+}
|
---|
13 |
|
---|
14 | interface
|
---|
15 |
|
---|
16 | uses
|
---|
17 | Classes, SysUtils, uos_dsp_utils ;
|
---|
18 |
|
---|
19 | type
|
---|
20 | PFFT = ^TFFT;
|
---|
21 |
|
---|
22 | { TFFT }
|
---|
23 |
|
---|
24 | TFFT = object
|
---|
25 | BitReversed: PInteger;
|
---|
26 | SinTable: PSingle;
|
---|
27 | Points: Integer;
|
---|
28 | FPCSinTable: array of Single;
|
---|
29 | function InitializeFFT(FFTLen: Integer): PFFT; static;
|
---|
30 | procedure EndFFT;
|
---|
31 | function GetFFT(FFTLen: Integer): PFFT; static;
|
---|
32 | procedure ReleaseFFT;
|
---|
33 | procedure InverseRealFFTf(buffer: PSingle);
|
---|
34 | procedure CleanupFFT; static; // ???
|
---|
35 | procedure RealFFTf(buffer: PSingle);
|
---|
36 | procedure ReorderToTime(Buffer: PSingle; TimeOut: PSingle);
|
---|
37 | procedure ReorderToFreq(Buffer: PSingle; RealOut: PSingle; ImagOut: PSingle);
|
---|
38 | end;
|
---|
39 |
|
---|
40 | type
|
---|
41 | PPSingle=^PSingle;
|
---|
42 | TSingleArray = array of Single;
|
---|
43 |
|
---|
44 | TNoiseRemoval = class;
|
---|
45 | TNoiseWriteProc = procedure(ASender: TObject; AData: PSingle; ASampleCount: Integer) of Object;
|
---|
46 |
|
---|
47 | { TNoiseRemoval }
|
---|
48 |
|
---|
49 | TNoiseRemoval = class
|
---|
50 | private
|
---|
51 | FDoProfile: Boolean;
|
---|
52 | FHasProfile: Boolean;
|
---|
53 |
|
---|
54 | // Parameters chosen before the first phase
|
---|
55 | FSampleRate: Double;
|
---|
56 | FWindowSize: Integer;
|
---|
57 | FSpectrumSize: Integer;
|
---|
58 | FMinSignalTime: Single; // in secs
|
---|
59 |
|
---|
60 | // The frequency-indexed noise threshold derived during the first
|
---|
61 | // phase of analysis
|
---|
62 | FNoiseThreshold: array of Single; // length in FSpectrumSize
|
---|
63 |
|
---|
64 | // Parameters that affect the noise removal, regardless of how the
|
---|
65 | // noise profile was extracted
|
---|
66 | FSensitivity: Double;
|
---|
67 | FFreqSmoothingHz: Double;
|
---|
68 | FNoiseGain: Double; // in dB, should be negative
|
---|
69 | FAttackDecayTime: Double;
|
---|
70 | FbLeaveNoise: Boolean;
|
---|
71 |
|
---|
72 | // change this later
|
---|
73 | procedure Initialize;
|
---|
74 | procedure Reset; // StartNewTrack
|
---|
75 | procedure ProcessSamples(len: Integer; Buffer:PSingle);
|
---|
76 | procedure FillFirstHistoryWindow;
|
---|
77 | procedure ApplyFreqSmoothing(ASpec: PSingle);
|
---|
78 | procedure GetProfile;
|
---|
79 | procedure RemoveNoise;
|
---|
80 | procedure RotateHistoryWindows;
|
---|
81 | procedure FinishTrack;
|
---|
82 | procedure Cleanup;
|
---|
83 | private
|
---|
84 | //FOutputTrack: PSingle; // WaveTrack;
|
---|
85 | FInSampleCount: Integer;
|
---|
86 | FOutSampleCount: Integer;
|
---|
87 | FInputPos: Integer;
|
---|
88 | FFFT: PFFT;
|
---|
89 | FFFTBuffer: PSingle; // FWindowSize
|
---|
90 | FWindow: PSingle; // FWindowSize
|
---|
91 | FFreqSmoothingBins: Integer;
|
---|
92 | FAttackDecayBlocks: Integer;
|
---|
93 | FOneBlockAttackDecay: Single;
|
---|
94 | FNoiseAttenFactor: Single;
|
---|
95 | FSensitivityFactor: Single;
|
---|
96 | FMinSignalBlocks: Integer;
|
---|
97 | FHistoryLen: Integer;
|
---|
98 | FInWaveBuffer: PSingle; // FWindowSize
|
---|
99 | FOutOverlapBuffer: PSingle; // FWindowSize
|
---|
100 | FSpectrums: array of PSingle; // FHistoryLen x FSpectrumSize
|
---|
101 | FGains: array of PSingle; // FHistoryLen x FSpectrumSize
|
---|
102 | FRealFFTs: array of PSingle; // FHistoryLen x FWindowSize
|
---|
103 | FImagFFTs: array of PSingle; // FHistoryLen x FWindowSize
|
---|
104 | FWriteProc: TNoiseWriteProc;
|
---|
105 | FInited: Boolean;
|
---|
106 | FTotalRead: QWord;
|
---|
107 | function GetNoiseProfile: TSingleArray;
|
---|
108 | procedure SetAttackDecayTime(AValue: Double);
|
---|
109 | procedure SetFreqSmoothingHz(AValue: Double);
|
---|
110 | procedure SetGain(AValue: Double);
|
---|
111 | procedure SetNoiseProfile(AValue: TSingleArray);
|
---|
112 | procedure SetSensitivity(AValue: Double);
|
---|
113 | public
|
---|
114 | constructor Create;
|
---|
115 | destructor Destroy; override;
|
---|
116 | function Init(ASampleRate: Integer): Boolean;
|
---|
117 | function Process(AData: PSingle; ASampleCount: Integer; AGetNoiseProfile: Boolean; AMoreComing: Boolean = False): Boolean;
|
---|
118 | procedure Flush; // finish writing out data in buffers.
|
---|
119 |
|
---|
120 | property NoiseProfile: TSingleArray read GetNoiseProfile write SetNoiseProfile;
|
---|
121 |
|
---|
122 | // these have defaults
|
---|
123 | property Sensitivity: Double read FSensitivity write SetSensitivity;
|
---|
124 | property Gain: Double read FNoiseGain write SetGain;
|
---|
125 | property FreqSmoothingHz: Double read FFreqSmoothingHz write SetFreqSmoothingHz;
|
---|
126 | property AttackDecayTime: Double read FAttackDecayTime write SetAttackDecayTime;
|
---|
127 | property LeaveNoise: Boolean read FbLeaveNoise write FbLeaveNoise;
|
---|
128 |
|
---|
129 | // don't mess with these.
|
---|
130 | property SampleRate: Double read FSampleRate;// write FSampleRate;
|
---|
131 | property WindowSize: Integer read FWindowSize;// write FWindowSize;
|
---|
132 | property SpectrumSize: Integer read FSpectrumSize;// write FSpectrumSize;
|
---|
133 | property MinSignalTime: Single read FMinSignalTime;// write FMinSignalTime; // in secs
|
---|
134 |
|
---|
135 | // This must be assigned or av's will occur
|
---|
136 | property WriteProc: TNoiseWriteProc read FWriteProc write FWriteProc;
|
---|
137 | end;
|
---|
138 |
|
---|
139 | type
|
---|
140 |
|
---|
141 | { TNoiseRemovalChannel }
|
---|
142 |
|
---|
143 | TNoiseRemovalChannel = class(TNoiseRemoval, IPAIODataIOInterface)
|
---|
144 | HasProfile: Boolean;
|
---|
145 | ProfileComplete: Boolean;
|
---|
146 | procedure WriteDataIO(ASender: IPAIODataIOInterface; AData: PSingle; ASamples: Integer);
|
---|
147 | end;
|
---|
148 |
|
---|
149 | { TNoiseRemovalMultiChannel }
|
---|
150 |
|
---|
151 | TNoiseRemovalMultiChannel = class(IPAIODataIOInterface)
|
---|
152 | private
|
---|
153 | FChannels,
|
---|
154 | FSampleRate: Integer;
|
---|
155 | FHelper: TPAIOChannelHelper;
|
---|
156 | FNoise: array of TNoiseRemovalChannel;
|
---|
157 | FWriteProc: TNoiseWriteProc;
|
---|
158 | //IPAIODataIOInterface
|
---|
159 | procedure WriteDataIO(ASender: IPAIODataIOInterface; AData: PSingle; ASamples: Integer);
|
---|
160 | procedure DataWrite(ASender: TObject; AData: PSingle; ASampleCount: Integer);
|
---|
161 | public
|
---|
162 | constructor Create(AChannels: Integer; ASampleRate: Integer);
|
---|
163 | destructor Destroy; override;
|
---|
164 | procedure ReadNoiseProfile(AData: PSingle; ASamples: Integer);
|
---|
165 | procedure ProcessNoise(AData: PSingle; ASamples: Integer);
|
---|
166 | procedure Flush;
|
---|
167 | property WriteProc: TNoiseWriteProc read FWriteProc write FWriteProc;
|
---|
168 | end;
|
---|
169 |
|
---|
170 | { TuosNoiseRemoval }
|
---|
171 | type
|
---|
172 | TuosNoiseRemoval = class(TNoiseRemovalMultiChannel)
|
---|
173 | OutStream: TStream;
|
---|
174 | public
|
---|
175 | isprofiled : boolean ;
|
---|
176 | samprate : integer ;
|
---|
177 | procedure WriteData(ASender: TObject; AData: PSingle; ASampleCount: Integer) ;
|
---|
178 | function FilterNoise(ANoisyAudio: PSingle; InFrames: Integer; out Samples: Integer): PSingle;
|
---|
179 | end;
|
---|
180 |
|
---|
181 | implementation
|
---|
182 | uses
|
---|
183 | math;
|
---|
184 |
|
---|
185 | const
|
---|
186 | PI = 3.14159265358979323846;
|
---|
187 | MAX_HFFT = 10;
|
---|
188 | var
|
---|
189 | FFTArray: array[0..MAX_HFFT-1] of PFFT;
|
---|
190 | FFTLockCount: array[0..MAX_HFFT-1] of Integer;
|
---|
191 |
|
---|
192 | { TMultiChannelNoiseRemoval }
|
---|
193 |
|
---|
194 | procedure TNoiseRemovalMultiChannel.WriteDataIO(ASender: IPAIODataIOInterface; AData: PSingle; ASamples: Integer);
|
---|
195 | begin
|
---|
196 | if Assigned(FWriteProc) then
|
---|
197 | FWriteProc(Self, AData, ASamples);
|
---|
198 | end;
|
---|
199 |
|
---|
200 | procedure TNoiseRemovalMultiChannel.DataWrite(ASender: TObject; AData: PSingle; ASampleCount: Integer);
|
---|
201 | begin
|
---|
202 | (FHelper as IPAIODataIOInterface).WriteDataIO(ASender as IPAIODataIOInterface, AData, ASampleCount);
|
---|
203 | end;
|
---|
204 |
|
---|
205 | constructor TNoiseRemovalMultiChannel.Create(AChannels: Integer;
|
---|
206 | ASampleRate: Integer);
|
---|
207 | var
|
---|
208 | i: Integer;
|
---|
209 | begin
|
---|
210 | FChannels:=AChannels;
|
---|
211 | FSampleRate:=ASampleRate;
|
---|
212 | FHelper := TPAIOChannelHelper.Create(Self);
|
---|
213 | SetLength(FNoise, AChannels);
|
---|
214 | for i := 0 to High(FNoise) do
|
---|
215 | begin
|
---|
216 | FNoise[i] := TNoiseRemovalChannel.Create;
|
---|
217 | FNoise[i].WriteProc:=@DataWrite;
|
---|
218 | FNoise[i].Init(ASampleRate);
|
---|
219 | FHelper.Outputs.Add(FNoise[i] as IPAIODataIOInterface);
|
---|
220 | end;
|
---|
221 | end;
|
---|
222 |
|
---|
223 | destructor TNoiseRemovalMultiChannel.Destroy;
|
---|
224 | var
|
---|
225 | i: Integer;
|
---|
226 | begin
|
---|
227 | for i := 0 to High(FNoise) do
|
---|
228 | begin
|
---|
229 | FNoise[i].Free;
|
---|
230 | end;
|
---|
231 | SetLength(FNoise, 0);
|
---|
232 | FHelper.Free;
|
---|
233 | end;
|
---|
234 |
|
---|
235 | procedure TNoiseRemovalMultiChannel.ReadNoiseProfile(AData: PSingle;
|
---|
236 | ASamples: Integer);
|
---|
237 | var
|
---|
238 | i: Integer;
|
---|
239 | begin
|
---|
240 | FHelper.Write(AData, ASamples);
|
---|
241 | for i := 0 to High(FNoise) do
|
---|
242 | begin
|
---|
243 | FNoise[i].ProfileComplete:=True;
|
---|
244 | FNoise[i].Process(nil, 0, True, False);
|
---|
245 | FNoise[i].HasProfile:=True;
|
---|
246 | FNoise[i].Init(FSampleRate);
|
---|
247 | end;
|
---|
248 | end;
|
---|
249 |
|
---|
250 | procedure TNoiseRemovalMultiChannel.ProcessNoise(AData: PSingle;
|
---|
251 | ASamples: Integer);
|
---|
252 | begin
|
---|
253 | FHelper.Write(AData, ASamples);
|
---|
254 | end;
|
---|
255 |
|
---|
256 | procedure TNoiseRemovalMultiChannel.Flush;
|
---|
257 | var
|
---|
258 | i: Integer;
|
---|
259 | begin
|
---|
260 | for i := 0 to High(FNoise) do
|
---|
261 | FNoise[i].Flush;
|
---|
262 | end;
|
---|
263 |
|
---|
264 | procedure TNoiseRemovalChannel.WriteDataIO(ASender: IPAIODataIOInterface;
|
---|
265 | AData: PSingle; ASamples: Integer);
|
---|
266 | begin
|
---|
267 | Process(AData, ASamples, not HasProfile, not HasProfile);
|
---|
268 | end;
|
---|
269 |
|
---|
270 | { TuosNoiseRemoval }
|
---|
271 |
|
---|
272 | procedure TuosNoiseRemoval.WriteData(ASender: TObject; AData: PSingle;
|
---|
273 | ASampleCount: Integer);
|
---|
274 | begin
|
---|
275 | OutStream.Write(AData^, ASampleCount*SizeOf(Single));
|
---|
276 | end;
|
---|
277 |
|
---|
278 | function TuosNoiseRemoval.FilterNoise(ANoisyAudio: PSingle; InFrames: Integer; out Samples: Integer): PSingle;
|
---|
279 | var
|
---|
280 | MNoisyAudio: TMemoryStream;
|
---|
281 | begin
|
---|
282 | OutStream := TMemoryStream.Create;
|
---|
283 |
|
---|
284 | MNoisyAudio := TMemoryStream.Create;
|
---|
285 | MNoisyAudio.Write(ANoisyAudio^, InFrames*SizeOf(Single));
|
---|
286 | MNoisyAudio.Position:=0;
|
---|
287 |
|
---|
288 | if isprofiled = false then // take the first chunk as noisy sample
|
---|
289 | begin
|
---|
290 | ReadNoiseProfile(PSingle(MNoisyAudio.Memory), MNoisyAudio.Size div SizeOf(Single));
|
---|
291 | isprofiled := true;
|
---|
292 | end;
|
---|
293 |
|
---|
294 | Result := nil;
|
---|
295 |
|
---|
296 | ProcessNoise(PSingle(MNoisyAudio.Memory), MNoisyAudio.Size div SizeOf(Single));
|
---|
297 |
|
---|
298 | Result:=GetMem(OutStream.Size);
|
---|
299 | Samples := OutStream.Size div SizeOf(Single);
|
---|
300 | OutStream.Position:=0;
|
---|
301 | OutStream.Read(Result^, OutStream.Size);
|
---|
302 |
|
---|
303 | MNoisyAudio.free;
|
---|
304 | OutStream.Free;
|
---|
305 |
|
---|
306 | // Result := ANoisyAudio;
|
---|
307 | end;
|
---|
308 |
|
---|
309 | { TNoiseRemoval }
|
---|
310 |
|
---|
311 | function NewFloatArray(ALength: Integer): PSingle; inline;
|
---|
312 | begin
|
---|
313 | Result := Getmem(ALength*SizeOf(Single));
|
---|
314 | end;
|
---|
315 |
|
---|
316 | procedure TNoiseRemoval.Initialize;
|
---|
317 | var
|
---|
318 | i: Integer;
|
---|
319 | begin
|
---|
320 | FFreqSmoothingBins := Trunc(FFreqSmoothingHz * FWindowSize / FSampleRate);
|
---|
321 | FAttackDecayBlocks := 1 + Trunc(FAttackDecayTime * FSampleRate / (FWindowSize / 2));
|
---|
322 | // Applies to amplitudes, divide by 20:
|
---|
323 | FNoiseAttenFactor := power(10, FNoiseGain/20);
|
---|
324 |
|
---|
325 | // Applies to gain factors which apply to amplitudes, divide by 20:
|
---|
326 | //FOneBlockAttackDecay := power(10.0, (FNoiseGain / (20.0 * FAttackDecayBlocks)));
|
---|
327 | FOneBlockAttackDecay := power(10.0, (FNoiseGain / FAttackDecayBlocks) / 20 );
|
---|
328 | // Applies to power, divide by 10:
|
---|
329 | FSensitivityFactor := power(10.0, FSensitivity/10.0);
|
---|
330 | FMinSignalBlocks := Trunc(FMinSignalTime * FSampleRate / (FWindowSize / 2));
|
---|
331 | if( FMinSignalBlocks < 1 ) then
|
---|
332 | FMinSignalBlocks := 1;
|
---|
333 | FHistoryLen := (2 * FAttackDecayBlocks) - 1;
|
---|
334 |
|
---|
335 | if (FHistoryLen < FMinSignalBlocks) then
|
---|
336 | FHistoryLen := FMinSignalBlocks;
|
---|
337 |
|
---|
338 | SetLength(FSpectrums, FHistoryLen);
|
---|
339 | SetLength(FGains, FHistoryLen);
|
---|
340 | SetLength(FRealFFTs, FHistoryLen);
|
---|
341 | SetLength(FImagFFTs, FHistoryLen);
|
---|
342 | for i := 0 to FHistoryLen-1 do
|
---|
343 | begin
|
---|
344 | FSpectrums[i] := NewFloatArray(FSpectrumSize);
|
---|
345 | FGains[i] := NewFloatArray(FSpectrumSize);
|
---|
346 | FRealFFTs[i] := NewFloatArray(FSpectrumSize);
|
---|
347 | FImagFFTs[i] := NewFloatArray(FSpectrumSize);
|
---|
348 | end;
|
---|
349 |
|
---|
350 | // Initialize the FFT
|
---|
351 | FFFT := TFFT.InitializeFFT(FWindowSize);
|
---|
352 |
|
---|
353 | FFFTBuffer := NewFloatArray(FWindowSize);
|
---|
354 | FInWaveBuffer := NewFloatArray(FWindowSize);
|
---|
355 | FWindow := NewFloatArray(FWindowSize);
|
---|
356 | FOutOverlapBuffer := NewFloatArray(FWindowSize);
|
---|
357 |
|
---|
358 | // Create a Hanning window function
|
---|
359 | for i := 0 to FWindowSize-1 do
|
---|
360 | FWindow[i] := 0.5 - 0.5 * cos((2.0*pi*i) / FWindowSize);
|
---|
361 |
|
---|
362 | if FDoProfile then
|
---|
363 | begin
|
---|
364 | FillChar(FNoiseThreshold[0], SizeOf(FNoiseThreshold[0])*FSpectrumSize, 0);
|
---|
365 | //for i := 0 to FSpectrumSize-1 do
|
---|
366 | // FNoiseThreshold[i] := float(0);
|
---|
367 | end;
|
---|
368 |
|
---|
369 | end;
|
---|
370 |
|
---|
371 | procedure TNoiseRemoval.Reset;
|
---|
372 | var
|
---|
373 | i, j: Integer;
|
---|
374 | begin
|
---|
375 | for i := 0 to FHistoryLen-1 do
|
---|
376 | begin
|
---|
377 | for j := 0 to FSpectrumSize-1 do
|
---|
378 | begin
|
---|
379 | FSpectrums[i][j] := 0;
|
---|
380 | FGains[i][j] := FNoiseAttenFactor;
|
---|
381 | FRealFFTs[i][j] := 0.0;
|
---|
382 | FImagFFTs[i][j] := 0.0;
|
---|
383 | end;
|
---|
384 | end;
|
---|
385 |
|
---|
386 | for j := 0 to FWindowSize-1 do
|
---|
387 | FOutOverlapBuffer[j] := 0.0;
|
---|
388 |
|
---|
389 | FInputPos := 0;
|
---|
390 | FInSampleCount := 0;
|
---|
391 | FOutSampleCount := -(FWindowSize div 2) * (FHistoryLen - 1);
|
---|
392 | end;
|
---|
393 |
|
---|
394 | function Min(A, B: Integer): Integer;
|
---|
395 | begin
|
---|
396 | if A < B then
|
---|
397 | Exit(A);
|
---|
398 | Result := B;
|
---|
399 | end;
|
---|
400 |
|
---|
401 | procedure TNoiseRemoval.ProcessSamples(len: Integer; Buffer: PSingle);
|
---|
402 | var
|
---|
403 | i: Integer;
|
---|
404 | avail: Integer;
|
---|
405 | begin
|
---|
406 | //while((len and FOutSampleCount) < FInSampleCount) do
|
---|
407 | while len > 0 do
|
---|
408 | begin
|
---|
409 | avail := Min(len, FWindowSize - FInputPos);
|
---|
410 | for i := 0 to avail-1 do
|
---|
411 | FInWaveBuffer[FInputPos + i] := buffer[i];
|
---|
412 | buffer += avail;
|
---|
413 | len -= avail;
|
---|
414 | FInputPos += avail;
|
---|
415 |
|
---|
416 | if (FInputPos = FWindowSize) then
|
---|
417 | begin
|
---|
418 | FillFirstHistoryWindow();
|
---|
419 | if (FDoProfile) then
|
---|
420 | GetProfile()
|
---|
421 | else
|
---|
422 | RemoveNoise();
|
---|
423 | RotateHistoryWindows();
|
---|
424 |
|
---|
425 | // Rotate halfway for overlap-add
|
---|
426 | //for(i = 0; i < mWindowSize / 2; i++) {
|
---|
427 | for i := 0 to FWindowSize div 2 -1 do
|
---|
428 | FInWaveBuffer[i] := FInWaveBuffer[i + FWindowSize div 2];
|
---|
429 |
|
---|
430 | FInputPos := FWindowSize div 2;
|
---|
431 | end;
|
---|
432 | end;
|
---|
433 | end;
|
---|
434 |
|
---|
435 | procedure TNoiseRemoval.FillFirstHistoryWindow;
|
---|
436 | var
|
---|
437 | i: Integer;
|
---|
438 | begin
|
---|
439 | for i := 0 to FWindowSize-1 do
|
---|
440 | FFFTBuffer[i] := FInWaveBuffer[i];
|
---|
441 | FFFT^.RealFFTf(FFFTBuffer);
|
---|
442 | //for(i = 1; i < (mSpectrumSize-1); i++) {
|
---|
443 | for i := 1 to FSpectrumSize-2 do
|
---|
444 | begin
|
---|
445 | FRealFFTs[0][i] := FFFTBuffer[FFFT^.BitReversed[i] ];
|
---|
446 | FImagFFTs[0][i] := FFFTBuffer[FFFT^.BitReversed[i]+1];
|
---|
447 | FSpectrums[0][i] := FRealFFTs[0][i]*FRealFFTs[0][i] + FImagFFTs[0][i]*FImagFFTs[0][i];
|
---|
448 | FGains[0][i] := FNoiseAttenFactor;
|
---|
449 | end;
|
---|
450 |
|
---|
451 | // DC and Fs/2 bins need to be handled specially
|
---|
452 | FSpectrums[0][0] := FFFTBuffer[0]*FFFTBuffer[0];
|
---|
453 | FSpectrums[0][FSpectrumSize-1] := FFFTBuffer[1]*FFFTBuffer[1];
|
---|
454 | FGains[0][0] := FNoiseAttenFactor;
|
---|
455 | FGains[0][FSpectrumSize-1] := FNoiseAttenFactor;
|
---|
456 | end;
|
---|
457 |
|
---|
458 | function Max(A,B: Integer): Integer; inline;
|
---|
459 | begin
|
---|
460 | if A>B then
|
---|
461 | Exit(A);
|
---|
462 | Result := B;
|
---|
463 | end;
|
---|
464 |
|
---|
465 | procedure TNoiseRemoval.ApplyFreqSmoothing(ASpec: PSingle);
|
---|
466 | var
|
---|
467 | tmp: PSingle;
|
---|
468 | i, j, j0, j1: Integer;
|
---|
469 | begin
|
---|
470 | tmp := NewFloatArray(FSpectrumSize);
|
---|
471 | for i := 0 to FSpectrumSize-1 do
|
---|
472 | begin
|
---|
473 | j0 := Max(0, i - FFreqSmoothingBins);
|
---|
474 | j1 := Min(FSpectrumSize-1, i + FFreqSmoothingBins);
|
---|
475 | tmp[i] := 0.0;
|
---|
476 | //for(j = j0; j <= j1; j++)
|
---|
477 | for j := j0 to j1-1 do
|
---|
478 | begin
|
---|
479 | tmp[i] += Aspec[j];
|
---|
480 | end;
|
---|
481 | tmp[i] := tmp[i] / (j1 - j0 + 1);
|
---|
482 | end;
|
---|
483 |
|
---|
484 | //for(i = 0; i < mSpectrumSize; i++)
|
---|
485 | for i := 0 to FSpectrumSize-1 do
|
---|
486 | Aspec[i] := tmp[i];
|
---|
487 |
|
---|
488 | Freemem(Tmp);
|
---|
489 | end;
|
---|
490 |
|
---|
491 | procedure TNoiseRemoval.GetProfile;
|
---|
492 | var
|
---|
493 | start,
|
---|
494 | finish,
|
---|
495 | i, j: Integer;
|
---|
496 | min: Single;
|
---|
497 | begin
|
---|
498 | // The noise threshold for each frequency is the maximum
|
---|
499 | // level achieved at that frequency for a minimum of
|
---|
500 | // mMinSignalBlocks blocks in a row - the max of a min.
|
---|
501 |
|
---|
502 | start := FHistoryLen - FMinSignalBlocks;
|
---|
503 | finish := FHistoryLen;
|
---|
504 |
|
---|
505 |
|
---|
506 | for j := 0 to FSpectrumSize-1 do
|
---|
507 | begin
|
---|
508 | min := FSpectrums[start][j];
|
---|
509 | for i := start+1 to finish-1 do
|
---|
510 | if (FSpectrums[i][j] < min) then
|
---|
511 | min := FSpectrums[i][j];
|
---|
512 |
|
---|
513 | if (min > FNoiseThreshold[j]) then
|
---|
514 | FNoiseThreshold[j] := min;
|
---|
515 | end;
|
---|
516 |
|
---|
517 | FOutSampleCount += FWindowSize div 2; // what is this for? Not used when we are getting the profile?
|
---|
518 | end;
|
---|
519 |
|
---|
520 | procedure TNoiseRemoval.RemoveNoise;
|
---|
521 | var
|
---|
522 | center: Integer;
|
---|
523 | start,
|
---|
524 | finish,
|
---|
525 | i,j : Integer;
|
---|
526 | min: Single;
|
---|
527 | out_: Integer;
|
---|
528 | begin
|
---|
529 | center := FHistoryLen div 2;
|
---|
530 | start := center - FMinSignalBlocks div 2;
|
---|
531 | finish := start + FMinSignalBlocks;
|
---|
532 |
|
---|
533 | // Raise the gain for elements in the center of the sliding history
|
---|
534 | for j := 0 to FSpectrumSize-1 do
|
---|
535 | begin
|
---|
536 | min := FSpectrums[start][j];
|
---|
537 | //for (i = start+1; i < finish; i++) {
|
---|
538 | for i := start+1 to finish-1 do
|
---|
539 | begin
|
---|
540 | if (FSpectrums[i][j] < min) then
|
---|
541 | min := FSpectrums[i][j];
|
---|
542 | end;
|
---|
543 | if (min > FSensitivityFactor * FNoiseThreshold[j]) and (FGains[center][j] < 1.0) then
|
---|
544 | begin
|
---|
545 | if (FbLeaveNoise) then
|
---|
546 | FGains[center][j] := 0.0
|
---|
547 | else
|
---|
548 | FGains[center][j] := 1.0;
|
---|
549 | end
|
---|
550 | else
|
---|
551 | begin
|
---|
552 | if (FbLeaveNoise) then
|
---|
553 | FGains[center][j] := 1.0;
|
---|
554 | end;
|
---|
555 | end;
|
---|
556 |
|
---|
557 | // Decay the gain in both directions;
|
---|
558 | // note that mOneBlockAttackDecay is less than 1.0
|
---|
559 | // of linear attenuation per block
|
---|
560 | for j := 0 to FSpectrumSize-1 do
|
---|
561 | begin
|
---|
562 | for i := center+1 to FHistoryLen-1 do
|
---|
563 | begin
|
---|
564 | if (FGains[i][j] < FGains[i - 1][j] * FOneBlockAttackDecay) then
|
---|
565 | FGains[i][j] := FGains[i - 1][j] * FOneBlockAttackDecay;
|
---|
566 | if (FGains[i][j] < FNoiseAttenFactor) then
|
---|
567 | FGains[i][j] := FNoiseAttenFactor;
|
---|
568 | end;
|
---|
569 | for i := center-1 downto 0 do
|
---|
570 | begin
|
---|
571 | if (FGains[i][j] < FGains[i + 1][j] * FOneBlockAttackDecay) then
|
---|
572 | FGains[i][j] := FGains[i + 1][j] * FOneBlockAttackDecay;
|
---|
573 | if (FGains[i][j] < FNoiseAttenFactor) then
|
---|
574 | FGains[i][j] := FNoiseAttenFactor;
|
---|
575 | end;
|
---|
576 | end;
|
---|
577 |
|
---|
578 |
|
---|
579 | // Apply frequency smoothing to output gain
|
---|
580 | out_ := FHistoryLen - 1; // end of the queue
|
---|
581 |
|
---|
582 | ApplyFreqSmoothing(FGains[out_]);
|
---|
583 |
|
---|
584 | // Apply gain to FFT
|
---|
585 | //for (j = 0; j < (mSpectrumSize-1); j++) {
|
---|
586 | for j := 0 to FSpectrumSize-2 do
|
---|
587 | begin
|
---|
588 | FFFTBuffer[j*2 ] := FRealFFTs[out_][j] * FGains[out_][j];
|
---|
589 | FFFTBuffer[j*2+1] := FImagFFTs[out_][j] * FGains[out_][j];
|
---|
590 | end;
|
---|
591 | // The Fs/2 component is stored as the imaginary part of the DC component
|
---|
592 | FFFTBuffer[1] := FRealFFTs[out_][FSpectrumSize-1] * FGains[out_][FSpectrumSize-1];
|
---|
593 |
|
---|
594 | // Invert the FFT into the output buffer
|
---|
595 | FFFT^.InverseRealFFTf(FFFTBuffer);
|
---|
596 |
|
---|
597 | // Overlap-add
|
---|
598 | for j := 0 to FSpectrumSize-2 do
|
---|
599 | begin
|
---|
600 | FOutOverlapBuffer[j*2 ] += FFFTBuffer[FFFT^.BitReversed[j] ] * FWindow[j*2 ];
|
---|
601 | FOutOverlapBuffer[j*2+1] += FFFTBuffer[FFFT^.BitReversed[j]+1] * FWindow[j*2+1];
|
---|
602 | end;
|
---|
603 |
|
---|
604 | // Output the first half of the overlap buffer, they're done -
|
---|
605 | // and then shift the next half over.
|
---|
606 |
|
---|
607 | if (FOutSampleCount >= 0) then // ...but not if it's the first half-window
|
---|
608 | begin
|
---|
609 | //FOutputTrack->Append((samplePtr)mOutOverlapBuffer, floatSample, mWindowSize / 2);
|
---|
610 | FWriteProc(Self, FOutOverlapBuffer, FWindowSize div 2);
|
---|
611 | end;
|
---|
612 |
|
---|
613 |
|
---|
614 | FOutSampleCount += FWindowSize div 2;
|
---|
615 | //for(j = 0; j < mWindowSize / 2; j++)
|
---|
616 | for j := 0 to FWindowSize div 2 -1 do
|
---|
617 | begin
|
---|
618 | FOutOverlapBuffer[j] := FOutOverlapBuffer[j + (FWindowSize div 2)];
|
---|
619 | FOutOverlapBuffer[j + (FWindowSize div 2)] := 0.0;
|
---|
620 | end
|
---|
621 | end;
|
---|
622 |
|
---|
623 | procedure TNoiseRemoval.RotateHistoryWindows;
|
---|
624 | var
|
---|
625 | last: Integer;
|
---|
626 | i: Integer;
|
---|
627 | lastSpectrum: PSingle;
|
---|
628 | lastGain: PSingle;
|
---|
629 | lastRealFFT: PSingle;
|
---|
630 | lastImagFFT: PSingle;
|
---|
631 | begin
|
---|
632 | last := FHistoryLen - 1;
|
---|
633 |
|
---|
634 | // Remember the last window so we can reuse it
|
---|
635 | lastSpectrum := FSpectrums[last];
|
---|
636 | lastGain := FGains[last];
|
---|
637 | lastRealFFT := FRealFFTs[last];
|
---|
638 | lastImagFFT := FImagFFTs[last];
|
---|
639 |
|
---|
640 | // Rotate each window forward
|
---|
641 | //for(i = last; i >= 1; i--) {
|
---|
642 | for i := last downto 1 do
|
---|
643 | begin
|
---|
644 | FSpectrums[i] := FSpectrums[i-1];
|
---|
645 | FGains[i] := FGains[i-1];
|
---|
646 | FRealFFTs[i] := FRealFFTs[i-1];
|
---|
647 | FImagFFTs[i] := FImagFFTs[i-1];
|
---|
648 | end;
|
---|
649 |
|
---|
650 | // Reuse the last buffers as the new first window
|
---|
651 | FSpectrums[0] := lastSpectrum;
|
---|
652 | FGains[0] := lastGain;
|
---|
653 | FRealFFTs[0] := lastRealFFT;
|
---|
654 | FImagFFTs[0] := lastImagFFT;
|
---|
655 | end;
|
---|
656 |
|
---|
657 | procedure TNoiseRemoval.FinishTrack;
|
---|
658 | var
|
---|
659 | empty: PSingle;
|
---|
660 | i: Integer;
|
---|
661 | begin
|
---|
662 | // Keep flushing empty input buffers through the history
|
---|
663 | // windows until we've output exactly as many samples as
|
---|
664 | // were input.
|
---|
665 | // Well, not exactly, but not more than mWindowSize/2 extra samples at the end.
|
---|
666 | // We'll delete them later in ProcessOne.
|
---|
667 | empty := NewFloatArray(FWindowSize div 2);
|
---|
668 | //for(i = 0; i < mWindowSize / 2; i++)
|
---|
669 | for i := 0 to FWindowSize div 2 -1 do
|
---|
670 | empty[i] := 0.0;
|
---|
671 |
|
---|
672 | while (FOutSampleCount < FInSampleCount) do
|
---|
673 | ProcessSamples(FWindowSize div 2, empty);
|
---|
674 |
|
---|
675 | Freemem(empty);
|
---|
676 | end;
|
---|
677 |
|
---|
678 | procedure TNoiseRemoval.Cleanup;
|
---|
679 | var
|
---|
680 | i: Integer;
|
---|
681 | begin
|
---|
682 | FFFT^.EndFFT;
|
---|
683 |
|
---|
684 | if (FDoProfile) then
|
---|
685 | ApplyFreqSmoothing(@FNoiseThreshold[0]);
|
---|
686 |
|
---|
687 |
|
---|
688 | for i := 0 to FHistoryLen-1 do
|
---|
689 | begin
|
---|
690 | FreeMem(FSpectrums[i]);
|
---|
691 | FreeMem(FGains[i]);
|
---|
692 | FreeMem(FRealFFTs[i]);
|
---|
693 | FreeMem(FImagFFTs[i]);
|
---|
694 | end;
|
---|
695 | SetLength(FSpectrums,0);
|
---|
696 | SetLength(FGains,0);
|
---|
697 | SetLength(FRealFFTs,0);
|
---|
698 | SetLength(FImagFFTs,0);
|
---|
699 |
|
---|
700 | FreeMem(FFFTBuffer);
|
---|
701 | FreeMem(FInWaveBuffer);
|
---|
702 | FreeMem(FWindow);
|
---|
703 | FreeMem(FOutOverlapBuffer);
|
---|
704 |
|
---|
705 | FInited := False;
|
---|
706 | end;
|
---|
707 |
|
---|
708 | function TNoiseRemoval.GetNoiseProfile: TSingleArray;
|
---|
709 | begin
|
---|
710 | SetLength(Result, FSpectrumSize);
|
---|
711 | Move(FNoiseThreshold[0], Result[0], FSpectrumSize);
|
---|
712 | end;
|
---|
713 |
|
---|
714 | procedure TNoiseRemoval.SetAttackDecayTime(AValue: Double);
|
---|
715 | begin
|
---|
716 | if FAttackDecayTime=AValue then Exit;
|
---|
717 | if AValue < 0.0 then AValue := 0;
|
---|
718 | if AValue > 1.0 then AValue := 1.0;
|
---|
719 | FAttackDecayTime:=AValue;
|
---|
720 | end;
|
---|
721 |
|
---|
722 | procedure TNoiseRemoval.SetFreqSmoothingHz(AValue: Double);
|
---|
723 | begin
|
---|
724 | if FFreqSmoothingHz=AValue then Exit;
|
---|
725 | if AValue<0 then AValue:=0;
|
---|
726 | if AValue>1000 then AValue := 1000;
|
---|
727 | FFreqSmoothingHz:=AValue;
|
---|
728 | end;
|
---|
729 |
|
---|
730 | procedure TNoiseRemoval.SetGain(AValue: Double);
|
---|
731 | begin
|
---|
732 | if FNoiseGain=AValue then Exit;
|
---|
733 | if AValue > 0 then AValue:=0;
|
---|
734 | if AValue < -48 then AValue := -48;
|
---|
735 |
|
---|
736 | FNoiseGain:=AValue;
|
---|
737 | end;
|
---|
738 |
|
---|
739 | procedure TNoiseRemoval.SetNoiseProfile(AValue: TSingleArray);
|
---|
740 | begin
|
---|
741 | SetLength(FNoiseThreshold, FSpectrumSize);
|
---|
742 | Move(AValue[0], FNoiseThreshold[0], FSpectrumSize);
|
---|
743 | FHasProfile:=True;
|
---|
744 |
|
---|
745 | FDoProfile:=False;
|
---|
746 | Cleanup; // set after FDoProfile so the profile is not processed.
|
---|
747 |
|
---|
748 | end;
|
---|
749 |
|
---|
750 | procedure TNoiseRemoval.SetSensitivity(AValue: Double);
|
---|
751 | begin
|
---|
752 | if FSensitivity=AValue then Exit;
|
---|
753 | if AValue < -20.0 then AValue:=-20.0;
|
---|
754 | if AValue > 20.0 then AValue := 20.0;
|
---|
755 | FSensitivity:=AValue;
|
---|
756 | end;
|
---|
757 |
|
---|
758 | constructor TNoiseRemoval.Create;
|
---|
759 | begin
|
---|
760 | FWindowSize:=2048;
|
---|
761 | FSpectrumSize:= 1 + FWindowSize div 2;
|
---|
762 |
|
---|
763 | // loaded prefs
|
---|
764 | FSensitivity := 0.0;
|
---|
765 | FNoiseGain := -24.0;
|
---|
766 | FFreqSmoothingHz := 150.0;
|
---|
767 | FAttackDecayTime:= 0.15;
|
---|
768 | FbLeaveNoise:=False;
|
---|
769 |
|
---|
770 | FMinSignalTime := 0.05;
|
---|
771 | FHasProfile := False;
|
---|
772 | FDoProfile := True;
|
---|
773 |
|
---|
774 | SetLength(FNoiseThreshold, FSpectrumSize);
|
---|
775 |
|
---|
776 | end;
|
---|
777 |
|
---|
778 | destructor TNoiseRemoval.Destroy;
|
---|
779 | begin
|
---|
780 | SetLength(FNoiseThreshold, 0);
|
---|
781 | inherited Destroy;
|
---|
782 | end;
|
---|
783 |
|
---|
784 | function TNoiseRemoval.Init(ASampleRate: Integer): Boolean;
|
---|
785 | begin
|
---|
786 | FSampleRate:=ASampleRate;
|
---|
787 | Initialize;
|
---|
788 | FInited:=True;
|
---|
789 | Result := True;
|
---|
790 | Reset;
|
---|
791 | end;
|
---|
792 |
|
---|
793 | function TNoiseRemoval.Process(AData: PSingle; ASampleCount: Integer;
|
---|
794 | AGetNoiseProfile: Boolean; AMoreComing: Boolean = False): Boolean;
|
---|
795 | begin
|
---|
796 | if not FInited then
|
---|
797 | Raise Exception.Create('TNoiseRemoval is not Inited');
|
---|
798 |
|
---|
799 | if not AGetNoiseProfile and not FHasProfile then
|
---|
800 | raise Exception.Create('Tried to remove noise without profile.');
|
---|
801 |
|
---|
802 | FDoProfile:=AGetNoiseProfile;
|
---|
803 |
|
---|
804 | if FDoProfile and (FTotalRead = 0) then
|
---|
805 | begin
|
---|
806 | Initialize;
|
---|
807 | reset;
|
---|
808 | end;
|
---|
809 |
|
---|
810 | Inc(FTotalRead, ASampleCount);
|
---|
811 | ProcessSamples(ASampleCount, AData);
|
---|
812 | Result := True;
|
---|
813 |
|
---|
814 | if AMoreComing then
|
---|
815 | Exit;
|
---|
816 |
|
---|
817 | if {AGetNoiseProfile or }FDoProfile then
|
---|
818 | begin
|
---|
819 | FHasProfile:=True;
|
---|
820 | Cleanup; // triggers the data in FNoiseThreshold to be processed
|
---|
821 |
|
---|
822 | // must be set after Cleanup() is called
|
---|
823 | //FDoProfile:=False;
|
---|
824 |
|
---|
825 | //Initialize;
|
---|
826 | FTotalRead:=0;
|
---|
827 | end;
|
---|
828 |
|
---|
829 | FHasProfile:=True;
|
---|
830 | FDoProfile := False;
|
---|
831 |
|
---|
832 | end;
|
---|
833 |
|
---|
834 | procedure TNoiseRemoval.Flush;
|
---|
835 | begin
|
---|
836 | if not FInited then
|
---|
837 | Exit;
|
---|
838 |
|
---|
839 | FinishTrack;
|
---|
840 | Cleanup; // Sets FInited to False
|
---|
841 | end;
|
---|
842 |
|
---|
843 | { TFFT }
|
---|
844 |
|
---|
845 | function TFFT.InitializeFFT(FFTLen: Integer): PFFT;
|
---|
846 | var
|
---|
847 | i: Integer;
|
---|
848 | temp: Integer;
|
---|
849 | mask: Integer;
|
---|
850 | begin
|
---|
851 | Result := New(PFFT);
|
---|
852 | if Result = nil then
|
---|
853 | raise EOutOfMemory.Create('Error allocating memory for FFT');
|
---|
854 |
|
---|
855 |
|
---|
856 | with Result^ do begin
|
---|
857 |
|
---|
858 | {*
|
---|
859 | * FFT size is only half the number of data points
|
---|
860 | * The full FFT output can be reconstructed from this FFT's output.
|
---|
861 | * (This optimization can be made since the data is real.)
|
---|
862 | *}
|
---|
863 | Points := FFTLen div 2;
|
---|
864 |
|
---|
865 | SetLength(FPCSinTable, 2*Points);
|
---|
866 | SinTable:=@FPCSinTable[0];
|
---|
867 |
|
---|
868 | BitReversed := Getmemory(Points*SizeOf(Integer));
|
---|
869 | if BitReversed = nil then
|
---|
870 | raise EOutOfMemory.Create('Error allocating memory for BitReversed.');
|
---|
871 |
|
---|
872 | for i := 0 to Points-1 do
|
---|
873 | begin
|
---|
874 | temp:=0;
|
---|
875 | mask := Points div 2;
|
---|
876 | while mask > 0 do
|
---|
877 | begin
|
---|
878 | //for(mask=h->Points/2;mask>0;mask >>= 1)
|
---|
879 | // temp=(temp >> 1) + (i&mask ? h->Points : 0);
|
---|
880 | temp := (temp shr 1);
|
---|
881 | if (i and mask) <> 0 then
|
---|
882 | temp := temp + Points;
|
---|
883 | //else temp := temp + 0; // why would you do that?
|
---|
884 | mask := mask shr 1;
|
---|
885 | end;
|
---|
886 |
|
---|
887 | BitReversed[i]:=temp;
|
---|
888 | end;
|
---|
889 |
|
---|
890 | for i := 0 to Points-1 do
|
---|
891 | begin
|
---|
892 | SinTable[BitReversed[i] ]:= -sin(2*PI*i/(2*Points));
|
---|
893 | SinTable[BitReversed[i]+1]:= -cos(2*PI*i/(2*Points));
|
---|
894 | end;
|
---|
895 |
|
---|
896 | {$ifdef EXPERIMENTAL_EQ_SSE_THREADED}
|
---|
897 | // new SSE FFT routines work on live data
|
---|
898 | for(i=0;i<32;i++)
|
---|
899 | if((1<<i)&fftlen)
|
---|
900 | h->pow2Bits=i;
|
---|
901 | {$endif}
|
---|
902 |
|
---|
903 | end; // with Result^
|
---|
904 |
|
---|
905 |
|
---|
906 | end;
|
---|
907 |
|
---|
908 | procedure TFFT.EndFFT;
|
---|
909 | begin
|
---|
910 | if Points>0 then
|
---|
911 | begin
|
---|
912 | Freemem(BitReversed);
|
---|
913 | SetLength(FPCSinTable, 0);
|
---|
914 | SinTable:=nil;
|
---|
915 | end;
|
---|
916 |
|
---|
917 | Dispose(PFFT(@Self));
|
---|
918 | end;
|
---|
919 |
|
---|
920 | function TFFT.GetFFT(FFTLen: Integer): PFFT;
|
---|
921 | var
|
---|
922 | h: Integer = 0;
|
---|
923 | n: Integer;
|
---|
924 | begin
|
---|
925 | n := fftlen div 2;
|
---|
926 |
|
---|
927 | while (h<MAX_HFFT) and (FFTArray[h] <> nil) and (n <> FFTArray[h]^.Points) do
|
---|
928 | begin
|
---|
929 | if (h<MAX_HFFT) then
|
---|
930 | begin
|
---|
931 | if(FFTArray[h] = nil) then
|
---|
932 | begin
|
---|
933 | FFTArray[h] := InitializeFFT(fftlen);
|
---|
934 | FFTLockCount[h] := 0;
|
---|
935 | end;
|
---|
936 | Inc(FFTLockCount[h]);
|
---|
937 | Exit(FFTArray[h]);
|
---|
938 | end
|
---|
939 | else begin
|
---|
940 | // All buffers used, so fall back to allocating a new set of tables
|
---|
941 | Exit(InitializeFFT(fftlen));
|
---|
942 | end;
|
---|
943 | Inc(h);
|
---|
944 | end;
|
---|
945 | end;
|
---|
946 |
|
---|
947 | procedure TFFT.ReleaseFFT;
|
---|
948 | var
|
---|
949 | h: Integer = 0;
|
---|
950 | begin
|
---|
951 |
|
---|
952 | while (h<MAX_HFFT) and (FFTArray[h] <> @Self) do
|
---|
953 | begin
|
---|
954 | if(h<MAX_HFFT) then
|
---|
955 | begin
|
---|
956 | Dec(FFTLockCount[h]);
|
---|
957 | end
|
---|
958 | else
|
---|
959 | begin
|
---|
960 | EndFFT;
|
---|
961 | end;
|
---|
962 | Inc(h);
|
---|
963 | end;
|
---|
964 | end;
|
---|
965 |
|
---|
966 | procedure TFFT.InverseRealFFTf(buffer: PSingle);
|
---|
967 | var
|
---|
968 | A, B: PSingle;
|
---|
969 | sptr: PSingle;
|
---|
970 | endptr1, endptr2: PSingle;
|
---|
971 | br1: PInteger;
|
---|
972 | HRplus,HRminus,HIplus,HIminus: Single;
|
---|
973 | v1,v2,sin,cos: Single;
|
---|
974 | ButterfliesPerGroup: Integer;
|
---|
975 | begin
|
---|
976 |
|
---|
977 | ButterfliesPerGroup:=Points div 2;
|
---|
978 |
|
---|
979 | //* Massage input to get the input for a real output sequence. */
|
---|
980 | A:=@buffer[2];
|
---|
981 | B:=@buffer[Points*2-2];
|
---|
982 | br1:=@BitReversed[1];
|
---|
983 | while(A<B) do
|
---|
984 | begin
|
---|
985 | sin:=SinTable[br1^];
|
---|
986 | cos:=SinTable[br1[1]];
|
---|
987 | //HRplus = (HRminus = *A - *B ) + (*B * 2);
|
---|
988 | HRminus:=A^-B^;
|
---|
989 | HRplus:=HRminus+ (B^ *2);
|
---|
990 |
|
---|
991 | //HIplus = (HIminus = *(A+1) - *(B+1)) + (*(B+1) * 2);
|
---|
992 | HIminus:=A[1]-B[1];
|
---|
993 | HIplus:=HIminus+(B[1] *2);
|
---|
994 |
|
---|
995 | v1 := (sin*HRminus + cos*HIplus);
|
---|
996 | v2 := (cos*HRminus - sin*HIplus);
|
---|
997 | A^ := (HRplus + v1) * single(0.5);
|
---|
998 | B^ := A^ - v1;
|
---|
999 | A[1] := (HIminus - v2) * single(0.5);
|
---|
1000 | B[1] := A[1] - HIminus;
|
---|
1001 |
|
---|
1002 | A+=2;
|
---|
1003 | B-=2;
|
---|
1004 | Inc(br1);
|
---|
1005 | end;
|
---|
1006 | //* Handle center bin (just need conjugate) */
|
---|
1007 | A[1] :=-A[1];
|
---|
1008 | {* Handle DC bin separately - this ignores any Fs/2 component
|
---|
1009 | buffer[1]=buffer[0]=buffer[0]/2;*}
|
---|
1010 | //* Handle DC and Fs/2 bins specially */
|
---|
1011 | //* The DC bin is passed in as the real part of the DC complex value */
|
---|
1012 | //* The Fs/2 bin is passed in as the imaginary part of the DC complex value */
|
---|
1013 | //* (v1+v2) = buffer[0] == the DC component */
|
---|
1014 | //* (v1-v2) = buffer[1] == the Fs/2 component */
|
---|
1015 | v1:=0.5*(buffer[0]+buffer[1]);
|
---|
1016 | v2:=0.5*(buffer[0]-buffer[1]);
|
---|
1017 | buffer[0]:=v1;
|
---|
1018 | buffer[1]:=v2;
|
---|
1019 |
|
---|
1020 | {*
|
---|
1021 | * Butterfly:
|
---|
1022 | * Ain-----Aout
|
---|
1023 | * \ /
|
---|
1024 | * / \
|
---|
1025 | * Bin-----Bout
|
---|
1026 | *}
|
---|
1027 |
|
---|
1028 | endptr1:=@buffer[Points*2];
|
---|
1029 |
|
---|
1030 | while(ButterfliesPerGroup>0) do
|
---|
1031 | begin
|
---|
1032 | A:=buffer;
|
---|
1033 | B:=@buffer[ButterfliesPerGroup*2];
|
---|
1034 | sptr:=@SinTable[0];
|
---|
1035 |
|
---|
1036 | while(A<endptr1) do
|
---|
1037 | begin
|
---|
1038 | sin:=sptr^; Inc(sptr); // *(sptr++);
|
---|
1039 | cos:=sptr^; Inc(sptr); // *(sptr++);
|
---|
1040 | endptr2:=B;
|
---|
1041 | while(A<endptr2) do
|
---|
1042 | begin
|
---|
1043 | v1:=B^*cos - B[1]*sin;
|
---|
1044 | v2:=B^*sin + B[1]*cos;
|
---|
1045 | B^ := (A^+v1)*Single(0.5);
|
---|
1046 | A^ := B^ - v1; Inc(A); Inc(B); //*(A++)=*(B++)-v1;
|
---|
1047 | B^ := (A^ + v2)* Single(0.5); //*B=(*A+v2)*(fft_type)0.5;
|
---|
1048 | A^ := B^ - v2; Inc(A); Inc(B); //*(A++)=*(B++)-v2;
|
---|
1049 | end;
|
---|
1050 | A:=B;
|
---|
1051 | B := @B[ButterfliesPerGroup*2];
|
---|
1052 | end;
|
---|
1053 | ButterfliesPerGroup := ButterfliesPerGroup shr 1;
|
---|
1054 | end;
|
---|
1055 | end;
|
---|
1056 |
|
---|
1057 | procedure TFFT.CleanupFFT;
|
---|
1058 | var
|
---|
1059 | h: Integer;
|
---|
1060 | begin
|
---|
1061 |
|
---|
1062 | for h :=0 to MAX_HFFT-1do begin
|
---|
1063 | if((FFTLockCount[h] <= 0) and (FFTArray[h] <> nil)) then
|
---|
1064 | begin
|
---|
1065 | FFTArray[h]^.EndFFT;
|
---|
1066 | FFTArray[h] := nil;
|
---|
1067 | end;
|
---|
1068 | end;
|
---|
1069 | end;
|
---|
1070 |
|
---|
1071 | procedure TFFT.RealFFTf(buffer: PSingle);
|
---|
1072 | var
|
---|
1073 | A, B: PSingle;
|
---|
1074 | sptr: PSingle;
|
---|
1075 | endptr1, endptr2: PSingle;
|
---|
1076 | br1, br2: PInteger;
|
---|
1077 | HRplus,HRminus,HIplus,HIminus: Single;
|
---|
1078 | v1,v2,sin_,cos_: Single;
|
---|
1079 | ButterfliesPerGroup: Integer;
|
---|
1080 | begin
|
---|
1081 | ButterfliesPerGroup:=Points div 2;
|
---|
1082 |
|
---|
1083 | {*
|
---|
1084 | * Butterfly:
|
---|
1085 | * Ain-----Aout
|
---|
1086 | * \ /
|
---|
1087 | * / \
|
---|
1088 | * Bin-----Bout
|
---|
1089 | *}
|
---|
1090 |
|
---|
1091 | endptr1:=buffer+Points*2;
|
---|
1092 |
|
---|
1093 | while(ButterfliesPerGroup>0) do
|
---|
1094 | begin
|
---|
1095 | A:=buffer;
|
---|
1096 | B:=buffer+ButterfliesPerGroup*2;
|
---|
1097 | sptr:=@SinTable[0];
|
---|
1098 |
|
---|
1099 | while(A<endptr1) do
|
---|
1100 | begin
|
---|
1101 | sin_:=sptr^;
|
---|
1102 | cos_ := sptr[1];
|
---|
1103 | endptr2:=B;
|
---|
1104 | while(A<endptr2) do
|
---|
1105 | begin
|
---|
1106 | v1 := B^ * cos_ + B[1] * sin_; //v1=*B*cos + *(B+1)*sin;
|
---|
1107 | v2 := B^ * sin_ - B[1] * cos_; //v2=*B*sin - *(B+1)*cos;
|
---|
1108 | B^ := A^+v1; //*B=(*A+v1);
|
---|
1109 | A^ := B^-2*v1; Inc(A); Inc(B); //*(A++)=*(B++)-2*v1;
|
---|
1110 |
|
---|
1111 | B^ := A^-v2; //*B=(*A-v2);
|
---|
1112 | A^ := B^+2*v2; Inc(A); Inc(B); //*(A++)=*(B++)+2*v2;
|
---|
1113 | end;
|
---|
1114 | A:=B;
|
---|
1115 | B:=B+ButterfliesPerGroup*2;
|
---|
1116 | sptr:=sptr+2;
|
---|
1117 | end;
|
---|
1118 |
|
---|
1119 | ButterfliesPerGroup := ButterfliesPerGroup shr 1;
|
---|
1120 | end;
|
---|
1121 |
|
---|
1122 | //* Massage output to get the output for a real input sequence. */
|
---|
1123 | br1:=@BitReversed[1]; // is this wrong? Should be @BitReversed[0] ; ?
|
---|
1124 | br2:=@BitReversed[Points-1];
|
---|
1125 |
|
---|
1126 | while(br1<br2) do
|
---|
1127 | begin
|
---|
1128 | sin_:=SinTable[br1[0]];
|
---|
1129 | cos_:=SinTable[br1[1]];
|
---|
1130 | A:=@buffer[br1^];
|
---|
1131 | B:=@buffer[br2^];
|
---|
1132 | //HRplus = (HRminus = *A - *B ) + (*B * 2);
|
---|
1133 | HRminus := A^ - B^;
|
---|
1134 | HRplus := HRminus + (B^ * 2);
|
---|
1135 |
|
---|
1136 | //HIplus = (HIminus = *(A+1) - *(B+1)) + (*(B+1) * 2);
|
---|
1137 | HIminus := A[1] - B[1];
|
---|
1138 | HIplus := HIminus + (B[1] * 2);
|
---|
1139 |
|
---|
1140 | v1 := (sin_*HRminus - cos_*HIplus);
|
---|
1141 | v2 := (cos_*HRminus + sin_*HIplus);
|
---|
1142 | A^ := (HRplus + v1) * single(0.5);
|
---|
1143 | B^ := A^ - v1;
|
---|
1144 | A[1] := (HIminus + v2) * single(0.5);
|
---|
1145 | B[1] := A[1] - HIminus;
|
---|
1146 |
|
---|
1147 | Inc(br1);
|
---|
1148 | Dec(br2);
|
---|
1149 | end;
|
---|
1150 | //* Handle the center bin (just need a conjugate) */
|
---|
1151 | A:=buffer+br1[1];
|
---|
1152 | A^ := -A^;
|
---|
1153 | {* Handle DC bin separately - and ignore the Fs/2 bin
|
---|
1154 | buffer[0]+=buffer[1];
|
---|
1155 | buffer[1]=(fft_type)0;*}
|
---|
1156 | ///* Handle DC and Fs/2 bins separately */
|
---|
1157 | ///* Put the Fs/2 value into the imaginary part of the DC bin */
|
---|
1158 | v1:=buffer[0]-buffer[1];
|
---|
1159 | buffer[0]+=buffer[1];
|
---|
1160 | buffer[1]:=v1;
|
---|
1161 |
|
---|
1162 | end;
|
---|
1163 |
|
---|
1164 | procedure TFFT.ReorderToTime(Buffer: PSingle; TimeOut: PSingle);
|
---|
1165 | var
|
---|
1166 | i: Integer;
|
---|
1167 | begin
|
---|
1168 | // Copy the data into the real outputs
|
---|
1169 | //for(int i=0;i<hFFT->Points;i++) {
|
---|
1170 | for i := 0 to Points-1 do
|
---|
1171 | begin
|
---|
1172 | TimeOut[i*2 ]:=buffer[BitReversed[i] ];
|
---|
1173 | TimeOut[i*2+1]:=buffer[BitReversed[i]+1];
|
---|
1174 | end;
|
---|
1175 | end;
|
---|
1176 |
|
---|
1177 | procedure TFFT.ReorderToFreq(Buffer: PSingle; RealOut: PSingle; ImagOut: PSingle
|
---|
1178 | );
|
---|
1179 | var
|
---|
1180 | i: Integer;
|
---|
1181 | begin
|
---|
1182 | // Copy the data into the real and imaginary outputs
|
---|
1183 | //for(int i=1;i<hFFT->Points;i++)
|
---|
1184 | for i := 1 to Points-1 do
|
---|
1185 | begin
|
---|
1186 |
|
---|
1187 | RealOut[i]:=buffer[BitReversed[i] ];
|
---|
1188 | ImagOut[i]:=buffer[BitReversed[i]+1];
|
---|
1189 | end;
|
---|
1190 | RealOut[0] := buffer[0]; // DC component
|
---|
1191 | ImagOut[0] := 0;
|
---|
1192 | RealOut[Points] := buffer[1]; // Fs/2 component
|
---|
1193 | ImagOut[Points] := 0;
|
---|
1194 | end;
|
---|
1195 |
|
---|
1196 | initialization
|
---|
1197 | FillChar(FFTArray, SizeOf(FFTArray), 0);
|
---|
1198 | FillChar(FFTLockCount, SizeOf(FFTLockCount), 0);
|
---|
1199 |
|
---|
1200 | end.
|
---|
1201 |
|
---|